Multivariate Analysis of Underwater Sounds

Submitted by

Kenneth John Powell

to the University of Exeter as a thesis for the degree of Doctor of Philosophy in Statistics in the Faculty of Science, May 1997

This thesis is available for Library use on the understanding that it is copyright material and that no quotation from the thesis may be published without proper acknowledgement.

I certify that all the material in this thesis which is not my own work has been identified and that no material is included for which a degree has previously been conferred upon me.

.....................

Abstract

This thesis considers the use of multivariate statistical methods in relation to a common signal processing problem, that of detecting features in sound recordings which contain interference, distortion and background noise.

Two separate but related areas of study are undertaken, first, the compression and noise reduction of sounds; second, the detection of intermittent departures (‘signals’) from the background sound environment (‘noise’), where the latter may be evolving and changing over time.

Compression and noise reduction are two closely related areas that have been studied for a wide range of signals, both one dimensional (such as sound) and two dimensional (such as images). Many well known techniques used in this field are based on the Fourier transform. In this work, we show how the comparative recent wavelet transform is superior for sound data involving short duration signals (such as shrimp clicks) whilst being at least as good as the Fourier transform for longer duration signals (such as dolphin whistles). Various noise reduction techniques involving thresholding wavelet transforms are examined and compared. We show how none of the standard threshold methods cope well with underwater sounds to any reasonable degree and propose a new technique, known as RunsThresh to overcome the perceived problems. The performance of this new method is contrasted with that of various standard thresholds.

Signal identification for underwater sounds is an area that has been examined in detail in much previous work. Here, we build upon the results gleaned from noise reduction to develop a methodology for detecting signals. The underwater noise environment is dynamically modelled using recursive density estimation of certain summary features of its wavelet decomposition. Observations which are considered to be outliers from this distribution are flagged as ‘signal’. The performance of our signal detection method is illustrated on artificial data, containing known signals, and on real data. This performance is compared with standard Fourier based methods for both cases.

Finally in this thesis, several ideas are presented and discussed which consider how the noise reduction and signal detection techniques examined in earlier chapters could be developed further, for example, in order to classify detected signals into different classes. These ideas are presented in outline only and are not followed up in detail, since they represent interesting directions for future study, rather than a primary focus of this thesis.

Acknowledgements

I would like to thank my supervisors, Prof. Wojtek Krzanowski and Dr. Trevor Bailey, for their help and support throughout this work. I would also like to acknowledge support for this research in the form of a Co-operative Award in Science and Engineering (CASE) studentship funded jointly by the Engineering and Physical Sciences Research Council (EPSRC) and DRA, Winfrith who also deserve a vote of thanks for supplying the data sets used in this research. I would especially like to thank Dr. Guy Nason of Bristol University who allowed me advanced access to WaveThresh 3 and his forthcoming book on wavelet analysis.

A general vote of thanks goes to all the staff and postgraduates in the department for making my time there enjoyable. Thanks also to those in Exeter who played sterling roles in keeping me sane, especially Aimée, Tom, Nigel, Fran and Hilary. Finally, I would like to thank my family for supporting me through this work.

Contents

1 Introduction 11

1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.1.1 Terminology . . . . . . . . . . . . . . . . . . . . . . . . 12 1.1.2 Introduction to the Data . . . . . . . . . . . . . . . . . 13

1.2 General Approaches . . . . . . . . . . . . . . . . . . . . . . . . 14 1.3 Structure of the Work . . . . . . . . . . . . . . . . . . . . . . 16

2 Review of Signal Processing 20 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.2 Signals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.2.1 Analogue Signals . . . . . . . . . . . . . . . . . . . . . 22 2.2.2 Digital Signals . . . . . . . . . . . . . . . . . . . . . . . 23

2.3 Signal Processing . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.3.1 Digital Filters . . . . . . . . . . . . . . . . . . . . . . . 28 2.3.2 The Fourier Transform . . . . . . . . . . . . . . . . . . 31 2.3.3 The Wavelet Transform . . . . . . . . . . . . . . . . . 42

2.4 Reasons for Using Wavelet Analysis . . . . . . . . . . . . . . . 56

3 About the Data 58 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.2 Technical Aspects of Underwater Sound Data . . . . . . . . . 59

3.2.1 Recording Circumstances and Data Format . . . . . . 59 3.2.2 Statistical Description of Underwater Sound . . . . . . 60

3.3 The Features of Underwater Sound . . . . . . . . . . . . . . . 60 3.3.1 Features Inherent in Underwater Sound . . . . . . . . . 61 3.3.2 Dolphin Signals . . . . . . . . . . . . . . . . . . . . . . 62 3.3.3 Shrimp Signals . . . . . . . . . . . . . . . . . . . . . . 66 3.3.4 Seal Signals . . . . . . . . . . . . . . . . . . . . . . . . 67 3.3.5 Falling Rain . . . . . . . . . . . . . . . . . . . . . . . . 67

3.4 The Features of the Artificial Sound . . . . . . . . . . . . . . . 72 3.4.1 Descriptive Analysis . . . . . . . . . . . . . . . . . . . 73

3.4.2 Usage of the Operator Data . . . . . . . . . . . . . . . 73

4 Data Compression And Noise Reduction 79 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.2 Data reduction using wavelets . . . . . . . . . . . . . . . . . . 80

4.2.1 Threshold Methods . . . . . . . . . . . . . . . . . . . . 81 4.2.2 RunsThresh . . . . . . . . . . . . . . . . . . . . . . . . 90

4.3 Comparisons . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

5 RunsThresh - Parameters and Effects 112 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 5.2 The RunsThresh Effect . . . . . . . . . . . . . . . . . . . . . . 112 5.3 The RunsThresh Parameters . . . . . . . . . . . . . . . . . . . 116

5.3.1 κ, the threshold effect . . . . . . . . . . . . . . . . . . 118 5.3.2 The Significance Level . . . . . . . . . . . . . . . . . . 121 5.3.3 The Block Size . . . . . . . . . . . . . . . . . . . . . . 123 5.3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 126

6 Signal Detection 127 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 6.2 Squaring the Pyramid . . . . . . . . . . . . . . . . . . . . . . 129

6.2.1 The Stationary Wavelet Transform . . . . . . . . . . . 129 6.2.2 Summarising the DWT . . . . . . . . . . . . . . . . . . 131

6.3 Signal Identification . . . . . . . . . . . . . . . . . . . . . . . . 134 6.3.1 Kernel Density Estimation . . . . . . . . . . . . . . . . 135 6.3.2 Identifying Sounds Using Probability Density Estimates 141 6.3.3 The Role of Thresholding in Identification . . . . . . . 150

6.4 Results and Comparisons . . . . . . . . . . . . . . . . . . . . . 151 6.4.1 The Fourier Technique . . . . . . . . . . . . . . . . . . 151

7 Conclusions and Further Work 167 7.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169

7.2.1 Classification . . . . . . . . . . . . . . . . . . . . . . . 170 7.2.2 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . 175

A Syntax for RunsThresh 176 A.1 Parameters for RunsThresh . . . . . . . . . . . . . . . . . . . 176 A.2 Example of Usage . . . . . . . . . . . . . . . . . . . . . . . . . 177 A.3 Credits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178 B Syntax for KP 179

B.1 Starting The Program . . . . . . . . . . . . . . . . . . . . . . 179 B.1.1 The Standard Matrix File Format . . . . . . . . . . . . 180 B.1.2 The Results File Format . . . . . . . . . . . . . . . . . 180

B.2 Parameters Used By KP . . . . . . . . . . . . . . . . . . . . . 180 B.3 Problems With Implementation . . . . . . . . . . . . . . . . . 182 B.4 Credits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183 B.5 The Program . . . . . . . . . . . . . . . . . . . . . . . . . . . 183

List of Figures

1.1 The time series for a dolphin whistle. . . . . . . . . . . . . . . 14 1.2 Frequency transforms of a dolphin whistle. . . . . . . . . . . . 17 1.3 The wavelet decomposition of a RunsThreshed dolphin whistle. 18 1.4 The reconstruction of a RunsThreshed dolphin whistle. . . . . 19

2.1 Various continuous time signals, after Ludeman (1986) …. 22 2.2 Various discrete time signals . . . . . . . . . . . . . . . . . . . 23 2.3 Sampling an analogue signal. . . . . . . . . . . . . . . . . . . . 24 2.4 Aliased signals. . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.5 Sampling positions. . . . . . . . . . . . . . . . . . . . . . . . . 27 2.6 Fourier transform of a signal. . . . . . . . . . . . . . . . . . . 28 2.7 Frequency responses of filters . . . . . . . . . . . . . . . . . . 30 2.8 Various window functions. . . . . . . . . . . . . . . . . . . . . 40 2.9 An example spectrogram. . . . . . . . . . . . . . . . . . . . . 41 2.10 Haar wavelet (ψ(t)) and the next detail level. . . . . . . . . . 45 2.11 Daubechies’ Extremal Phase order 4 wavelet. . . . . . . . . . . 48 2.12 Daubechies’ Extremal Phase order 10 wavelet. . . . . . . . . . 49 2.13 The Discrete Wavelet Transform for 16 data points. . . . . . . 51

3.1 Frequency decompositions of various signals. . . . . . . . . . . 63 3.2 Frequency decompositions of dolphin sounds. . . . . . . . . . . 64 3.3 Raw plots of dolphin sounds. . . . . . . . . . . . . . . . . . . . 65 3.4 Various representations of shrimp snaps. . . . . . . . . . . . . 68 3.5 Comparison of a shrimp snap and a dolphin click. . . . . . . . 69 3.6 Various representations of a seal bark. . . . . . . . . . . . . . 70 3.7 Various representations of rain falling. . . . . . . . . . . . . . 71 3.8 Various representations of the operator data at 5 dB NL. … 74 3.9 Operator data at 0 and −5 dB NL. . . . . . . . . . . . . . . . 75 3.10 Operator data at −8 and −10 dB SNR. . . . . . . . . . . . . . 76 3.11 The position of the word operator in the signal. . . . . . . . . 77

4.1 A dolphin whistle of length 0.8s. . . . . . . . . . . . . . . . . . 84 4.2 Wavelet decomposition of the whistle. . . . . . . . . . . . . . . 85 4.3 VisuShrink applied to the whistle. . . . . . . . . . . . . . . . . 86 4.4 SureShrink applied to the whistle. . . . . . . . . . . . . . . . . 87 4.5 Cross Validation applied to the whistle. . . . . . . . . . . . . . 89 4.6 Magnified section of the whistle wavelet decomposition. . . . . 90 4.7 RunsThresh applied to the whistle. . . . . . . . . . . . . . . . 92 4.8 The four test functions defined in Table 4.1 . . . . . . . . . . . 94 4.9 The four typical sounds used for the comparisons . . . . . . . 95 4.10 The four test functions with added noise . . . . . . . . . . . . 96 4.11 Test signals after RunsThresh has been applied . . . . . . . . 100 4.12 Test signals after VisuShrink has been applied . . . . . . . . . 101 4.13 Test signals after SureShrink has been applied . . . . . . . . . 102 4.14 Test signals after Cross-validation has been applied . . . . . . 103 4.15 Underwater signals after RunsThresh has been applied …. 106 4.16 Underwater signals after VisuShrink has been applied . . . . . 107 4.17 Underwater signals after SureShrink has been applied . . . . . 108 4.18 Underwater signals after Cross-validation has been applied . . 109

5.1 Significant blocks of coefficients (shown in black) during Run-

sThresh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.2 The wavelet decomposition of the whistle after thresholds have

been applied. . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

6.1 The process of sound identification. . . . . . . . . . . . . . . . 129 6.2 Two wavelet decompositions of a Doppler signal. . . . . . . . . 130 6.3 Kernel and data shapes if the Fukunaga method is used, after

Fukunaga (1972). . . . . . . . . . . . . . . . . . . . . . . . . . 136 6.4 The total energy by variable for various signals. . . . . . . . . 144 6.5 The total energy by a subset of variables for various signals. . 145 6.6 Results for the analysis of the rain data. . . . . . . . . . . . . 157 6.7 Results for the analysis of the seal data. . . . . . . . . . . . . 158 6.8 Results for the analysis of the shrimp data. . . . . . . . . . . . 159 6.9 Results for the analysis of the dolphin data. . . . . . . . . . . 160 6.10 The rain data before and after signal identification. . . . . . . 165 6.11 The seal data before and after signal identification. . . . . . . 165 6.12 The shrimp data before and after signal identification. . . . . . 166

B.1 The standard matrix file format . . . . . . . . . . . . . . . . . 180

List of Tables

4.1 Formulae for test functions . . . . . . . . . . . . . . . . . . . . 93 4.2 AMSE for standard signals using order 8 Least Asymmetric

wavelets with soft thresholding . . . . . . . . . . . . . . . . . . 97 4.3 AMSE for standard signals using order 2 Extremal Phase

wavelets with soft thresholding . . . . . . . . . . . . . . . . . . 98 4.4 AMSE for standard signals using order 8 Least Asymmetric

wavelets with hard thresholding . . . . . . . . . . . . . . . . . 99 4.5 AMSE for standard signals using order 2 Extremal Phase

wavelets with hard thresholding . . . . . . . . . . . . . . . . . 99 4.6 Quality of underwater sounds after thresholding using order 8

Least Asymmetric wavelets with soft thresholding . . . . . . . 104 4.7 Quality of underwater sounds after thresholding using order 8

Least Asymmetric wavelets with hard thresholding . . . . . . 104 4.8 Quality of underwater sounds after thresholding using order 2

Extremal Phase wavelets with soft thresholding . . . . . . . . 105 4.9 Quality of underwater sounds after thresholding using order 2

Extremal Phase wavelets with hard thresholding . . . . . . . . 105 4.10 Comparison of RunsThresh results using different wavelets . . 106 4.11 Compression ratios different thresholds methods . . . . . . . . 110 4.12 Compression ratios for different underwater sounds . . . . . . 111

5.1 Summary measures for dolphin whistles after thresholding us-

ing order 8 Least Asymmetric wavelets with soft thresholding 118 5.2 Summary measures for shrimp clicks after thresholding using

order 8 Least Asymmetric wavelets with soft thresholding … 119 5.3 Quality of dolphin whistles after thresholding using order 8

Least Asymmetric wavelets with soft thresholding . . . . . . . 120 5.4 Quality of shrimp clicks after thresholding using order 8 Least

Asymmetric wavelets with soft thresholding . . . . . . . . . . 120 5.5 Summary measures for dolphin whistles after thresholding us-

ing order 8 Least Asymmetric wavelets with soft thresholding 121 5.6 Summary measures for shrimp clicks after thresholding using

order 8 Least Asymmetric wavelets with soft thresholding … 122 5.7 Quality of dolphin whistles after thresholding using order 8

Least Asymmetric wavelets with soft thresholding . . . . . . . 123 5.8 Quality of shrimp clicks after thresholding using order 8 Least

Asymmetric wavelets with soft thresholding . . . . . . . . . . 123 5.9 Summary measures for dolphin whistles after thresholding us-

ing order 8 Least Asymmetric wavelets with soft thresholding 124 5.10 Summary measures for shrimp clicks after thresholding using

order 8 Least Asymmetric wavelets with soft thresholding … 125 5.11 Quality of dolphin whistles after thresholding using order 8

Least Asymmetric wavelets with soft thresholding . . . . . . . 125 5.12 Quality of shrimp clicks after thresholding using order 8 Least

Asymmetric wavelets with soft thresholding . . . . . . . . . . 126

6.1 Sample size required to produce acceptable mean square errors

in the estimation of the density at 0, after Silverman (1986). . 146 6.2 Error rates for the operator data using f̂(x_{n+1}) . . . . . . . . 154 6.3 Error rates for the operator data using δ(x_{n+1}, \overline{x}) . . . . . . . 154 6.4 Error rates for the operator data using δ(x_{n+1}, \overline{x_{w}}) . . . . . . . 155 6.5 Error rates for the operator data using \overline{δ(x_{n+1}, x_{t})} . . . . . . . 155 6.6 Error rates for the operator data using \frac{∑ f̂(x_{t})δ(x_{n+1},x_{t})}{∑ f̂(x_{t})} . . . . . 156 6.7 Error rates for the operator data using δ(x_{n+1}, x_{m}) . . . . . . 156 6.8 Error rates for the operator data using Fourier based energy

detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 6.9 Signal detection rates . . . . . . . . . . . . . . . . . . . . . . . 162 6.10 Comparisons of classification for rain data . . . . . . . . . . . 162 6.11 Comparisons of classification for seal data . . . . . . . . . . . 162 6.12 Comparisons of classification for shrimp data . . . . . . . . . . 163 6.13 Comparisons of classification for dolphin data . . . . . . . . . 163

Chapter 1

Introduction

1.1 Background

The issues examined in this work arose from practical problems experienced by the Defence Research Agency (DRA), Sonar and Underwater Technology Division, but there are applications in far wider areas such as geology, marine biology and structural engineering. At the most basic level, the question posed was “how can a computer examine real time sounds so that anomalous events, potential ‘signals’, are detected?”

In the context presented, for the purpose of this thesis, by the DRA, anomalous events are underwater animal signals whilst for structural engineering, for instance, such signals might be stress fractures forming in load bearing struts, or the movement of foundations before an earthquake.

The DRA’s motivation for solving this problem is the desire to acquire a system that can be used as a backup to the professional SONAR (SOund NAvigation and Ranging) personnel on board ships and submarines. Space is at a premium on the latter especially, and a system that means fewer staff are required would be most beneficial.

This work will proceed to examine the question posed and various possible solutions to it. The main body of the work will examine the problem by way of frequency^{1} transforms and non-parametric statistical methods. In advance of that, an overview of the specific underwater sound problem will be given along with a brief description of some of the methods previously used that are not considered in this work. This chapter will then finish with a summary of the remaining areas covered in this work. First however, the terminology used herein will be described.

1The frequency is the pitch of the sound.

1.1.1 Terminology

It is appropriate, before starting any explanations that some definitions of the more common terminology should be given. In this work, several common words are used in a very specific sense and particular mathematical notation is used in several places.

Any mention of sound, signal and noise throughout this thesis has a very specific meaning. A sound (sometimes referred to as raw data or original data) is the data collected by the hydrophone (underwater microphone) at source. In the context of this problem it is digital data collected at a sample rate measured in terms of cycles per second (Hz). For example, a sample rate of 8 kHz corresponds to 8000 observations per second. A signal is any part of the sound that is of interest. This could be, for instance, a dolphin whistle or a shrimp click. Signals can also be known as transients. Further examples of signals are given in chapter 3. Finally, noise is any part of the sound that is not of interest. This includes traditional background “noise” such as the sound made by the sea and also sounds produced, as part of the recording process or by artifacts, by objects such as motor boats that are not of interest for this problem.

From these definitions we can observe that a sound is actually a time series which can be split into two time series corresponding to the signal and the noise — the signal component can of course be absent at any particular time, and the noise component can change over time.

A general time series of length T is written as x(T) = x_{1}, x_{2}, …, x_{T −1}, x_{T}. This of course allows the signal, noise or sound to be written using the above notation. However, we have already seen that the sound is a combination of the noise and signal, hence sound(t) = signal_{1} + noise_{1}, signal_{2} + noise_{2}, …, signal_{T} + noise_{T}. Therefore it is easy to see that the noise is actually a combination of many time series (from the copious sources of noise in an underwater environment). Similarly the signal could potentially be made up from many sources both similar (multiple shrimps clicking) and different (dolphins barking and whistling). It should be noted that the combination of noise and signal is not always additive as implied here. Such combinations could, in many cases, be for instance, multiplicative.

1.1.2 Introduction to the Data

The data will be fully described in Chapter 3 so, at this point, it is sufficient to say that the underwater sound recordings used are of a similar quality to compact discs. The recordings were made around the shores of the British Isles, dolphins being recorded off South Wales and the shrimps recorded in the North Sea.

Figure 1.1 shows the raw time series for a dolphin whistle. This is a particularly clear example, hence the low level of noise running through the time period (the noise is the narrow central band, seen most clearly at the very beginning of the signal). With the data recorded at a rate of 40960 samples per second the data sets are clearly massive and constitute an unusually large amount of data to analyse. This is the first unusual feature of the data. The second such feature is that some signals studied are broad-band, that is, there is the potential to have signals including frequencies across the range from 0 kHz (also known as DC) through to at least 16 kHz.

Other unusual features of the data are the lack of identified test data sets (most signal identification techniques rely on databases of known signals with which to compare new sounds) and the existence of very short time signals or transients (measured in hundredths of a second).

These various features may be contrasted with the situation found in human speech analysis, which is the most common form of data used in sound identification and classification. Human speech has a much reduced frequency range. It is narrow band with a total range of the order of 4 kHz and a bias towards the lower end of the frequency spectrum. Also with speech data,

Figure 1.1: The time series for a dolphin whistle.

and indeed most signal identification, it is possible to build up a database of noise free, pre-classified example signals whereas with the underwater sounds analysed here this is not possible —dolphins do not respond well when asked to say “aaah”! Similarly, most forms of speech are relatively long whereas a shrimp or dolphin click is less than six hundredths of a second in length.

1.2 General Approaches

The primary goals of this work are noise reduction and signal identification. It should be noted in passing however that signal detection is often a necessary prerequisite to signal classification — discriminating between several different types of signal. This work therefore may have relevance to studies concerned with the classification of sonar and high-dimensional signals (such as images); either where neural networks are employed (e.g. Gorman and Sejnowski 1988; Ghosh et al. 1992; Ripley 1994; Smith et al. 1993); or where wavelets are used as a feature extraction method for discrimination or classification (e.g Coifman and Saito 1994; Telfer et al. 1994; Learned and Wilsky 1995). Some possible approaches to the classification problem will be discussed briefly as ideas for future work in the final chapter.

Clearly, there are many possible approaches to problems in the field of noise reduction and signal detection. Possible approaches may arise from such methods as filtering (Prasad et al. 1995), multivariate process control (Wierda 1994; Fuchs and Benjamini 1994), the detection of outliers in multivariate time series (Khattree and Naik 1987), multivariate discrimination techniques such as regression trees or principal component analysis (Krzanowski 1988; Atlas et al. 1990) or methods arising from the analysis of images or seismic signals (Chen 1985; Pal and Pal 1993).

Many of these methods for analysing sound data make use of the Fourier transform to allow the analysis to be performed in the frequency domain (Cohen 1989; Boashash 1990; Jones and Parks 1992). In this work, we also undertake the analysis in the frequency domain. We do not, however use the Fourier transform to undertake the transform from the time to the frequency domain. As will be shown in Chapter 2, the Fourier transform does not perform well in certain of the applications of interest to us, especially where very short time transients are involved. In this work, the wavelet transform will be used to perform the domain change. Although the use of wavelets for sound classification and compression is not new — for example, (Telfer et al. 1994) and (Learned and Wilsky 1995) use techniques which select or adapt wavelets from the wide range of basis that exist in the literature (Coifman and Saito 1994) — some of the techniques used in this work are especially adapted to cope with the unusual features of underwater sound.

Although wavelets are used in place of the Fourier transform, this still leaves the method of noise reduction or signal detection, in the frequency domain, to be decided. Traditionally, especially in the more engineering led approaches, some form of simple filtering, such as a Kalman filter (Mc-Causland 1969), would be used for the identification. In this work, wavelet thresholding methods are considered for noise reduction and a statistical approach for identification is developed using multivariate density estimation implemented via a multivariate kernel approach (Silverman 1986).

1.3 Structure of the Work

Given the background outlined in the previous two sections, this chapter will conclude with an overview of the remainder of this work. In detailing the contents of the remaining chapters, use will be made of the whistle shown in Figure 1.1 to demonstrate the different techniques used to analyse the sounds.

Chapter 2 gives a general introduction to digital signal processing. Definitions of digital and analogue signals will be given along with an introduction to filters and their uses in signal identification. Finally the two main timefrequency transforms, the Fourier and wavelet transforms, will be described.

The major topics introduced here are the Fourier transform and its associated representation, the spectrogram, and the wavelet transform with its associated wavelet decomposition plot. Figure 1.2 shows the dolphin whistle depicted in Figure 1.1 subjected to these two transforms.

Both graphs show time as the horizontal axis and frequency (or equivalent) as the vertical. The spectrogram shows high frequency at the top of the graph whilst the wavelet decomposition shows low frequency at the top. Full details of how these graphs are constructed are given in Chapter 2.

After this introduction, Chapter 3 examines the data sets used in detail. Descriptions are given of the format, the recording circumstances and the contents for the underwater sounds — recordings of dolphins, shrimps and rain. Finally a new data set consisting of human speech is introduced. This will be used for calibration and comparisons between methods.

Chapter 4 examines the different wavelet based techniques available for data reduction and denoising. Any real life data, be it from underwater or above ground, is likely to be contaminated with background noise. Audio playback would be greatly enhanced if such noise could be removed prior to presentation to a human listener. Indeed this is a useful first step in improving detection rates. There is also a good probability that there will be a requirement to store the incoming sounds for later analysis, with the data recorded at compact disc quality. This will obviously require large amounts of storage space per hour. It would be useful to be able to compress the stored

(a) Spectrogram of the whistle (b) Wavelet decomposition of the whistle

Figure 1.2: Frequency transforms of a dolphin whistle. data so that space requirements could be minimised, even if only certain time periods are preserved.

Chapter 4, then, examines various ways these two objectives can be achieved using the wavelet representation as a starting point. It introduces a modified thresholding technique to use on the transformed data which simultaneously denoises and compresses the data. This new technique is called RunsThresh and is described in detail. Figure 1.3 shows the wavelet decomposition of the whistle after RunsThreshing. In this case, a reduction of 39% has been achieved between the original representation and the thresholded version. Figure 1.4 shows the reconstruction of the whistle from the thresholded coefficients.

Figure 1.3: The wavelet decomposition of a RunsThreshed dolphin whistle.

More details of the RunsThresh method are given in Chapter 5. The different parameters that control the thresholding are examined and evidence is offered that RunsThresh is a significantly different technique.

The major problem of signal identification is dealt with in Chapter 6. This chapter explains how the wavelet decomposition is pre-processed before being passed into a multivariate kernel estimator. The probability density estimate

Figure 1.4: The reconstruction of a RunsThreshed dolphin whistle.

obtained is then used as a basis for describing the background noise. New observations are compared with this to determine whether they are noise or signal. This approach uses various methods to decide whether or not new observations are from the density. These methods are described in some detail, along with other possible approaches to this problem which have been developed elsewhere.

Finally, Chapter 7 summarises the conclusions gleaned from this work and attempts to show how the noise reduction and signal identification can be combined to produce a “clean” output for audio playback. It also discusses possible extensions to the approaches contained herein that would allow events to be classified by type as well as simply identified.

Review of Signal Processing

2.1 Introduction

Signal processing is a well established field whose practitioners are mainly in the field of engineering but which is also an area of major interest in statistics. At the most basic level, it is concerned with operating on a sound (or image) to extract useful information, known as the signal (Smith and Mersereau 1992). The operation may be a non-destructive transformation of the given data, but some important processing methods are destructive and so cannot be reversed.

There are various reasons for performing this operation. Features of the data may be easier to observe in the transformed space, for example some features of speech are more easily observed in the frequency rather than in the time domain; edges in an image can be visually enhanced by means of a filter, and filtering can give storage benefits in that the transformed sound can be stored more compactly. With music and speech, where the signal is known to be in certain frequency bands, all other frequencies can be filtered out, leaving a much clearer signal to work with.

Overall then, in most sound analysis problems, one of the first steps is to filter the data in such a way that the required features of the signal are enhanced. This chapter will discuss current signal processing practice in both the statistical and engineering fields (although there is much overlap between the two). A brief discussion of signal types will be given first, followed by a review of various methods of signal processing, with particular emphasis on digital filters, the Fourier transform and the wavelet transform. The chapter will end with a justification of the use of the wavelet transform in the remainder of this work.

2.2 Signals

Signals can come in various forms. There are one dimensional signals (often speech, SONAR, RADAR (RAdio Detecting And Ranging) or radio waves) and two dimensional signals (often images) both of which can be analogue (continuous time), digital (discrete time) or a combination of the two (such as sampled analogue data, as on a compact disc).

The underwater sound data analysed in this work is sampled analogue data, which has the properties of digital data — including some sampling problems discussed below. The main emphasis of this chapter is therefore on digital signal processing but a brief discussion of analogue signals will be given first as the data were originally of this form.

It is, within limits, always possible to convert from an analogue to a digital signal and back by means of A/D and D/A converters. The most basic form of an A/D (analogue to digital) converter is a simple sampler, whilst a D/A (digital to analogue) converter is a continuous output constructor and smoother in one unit. Details of these two units are not really necessary for what follows, but it is important to know they exist as they allow all sounds to be processed digitally. This has many advantages, some of which are detailed in the section on signal processing. Examples of an A/D converter are a sampler such as those used in many computers and a DAT (Digital Audio Tape) recorder. D/A devices include a compact disc player and computer sound cards.

2.2.1 Analogue Signals

An analogue signal is a continuous time function, which is assigned a real value for every value of time. Some of the more common analogue signals are shown in figure 2.1. The periodic signal might look to be discrete as it only takes two values and has large discontinuities but it is, in fact, analogue as it takes a real value (although from only two possibilities) for every value of time.

(a) Periodic  (b) Positive time transient

(c) Sinusoidal  (d) Random

Figure 2.1: Various continuous time signals, after Ludeman (1986)

The example of the random signal is similar to the original underwater signals used throughout this work. It is, obviously, the hardest form of signal to model.

In general, we will refer to a continuous time signal as x(t) to distinguish it from discrete time signals (which will be indexed by n).

2.2.2 Digital Signals

Digital signals are discrete time functions, that is, they are defined on the positive and negative integers. A discrete time signal is referred to as x_{n} for n = −∞, …, −1, 0, 1, …, ∞. Such a signal is periodic if an integer P > 0 exists such that x_{n+P} = x_{n} for all n. For signals that satisfy this condition, P is the period of the signal.

For a discrete time signal to be a digital signal, the values of x_{n} must take on a finite number of different values; if they take on a countably infinite set of values then they are not digital. Examples of discrete time and digital signals are shown in figure 2.2. Only w_{n} and z_{n} in this diagram are digital.

(a) Periodic digital  (b) Positive n

(c) Discrete sinusoidal  (d) Non-periodic digital

Figure 2.2: Various discrete time signals

The underwater signals dealt with throughout this work are digital as they take a value from the set {−32768, −32767, …, −1, 0, 1, …, 32766, 32767} every 1/40960 seconds.

Sampling

As mentioned above, analogue signals can be converted into their digital counterparts by means of an A/D converter. The main part of such a device is the sampler. Figure 2.3 shows a typical sampling process (Grant et al. 1989). An analogue signal is shown at the top, with a sampling waveform, a uniform sequence of equally spaced impulses, below. The sampled signal is then the product of these two signals. This sampled signal could be reconstructed with a D/A converter to retrieve the original waveform, or a close approximation to it.

Figure 2.3: Sampling an analogue signal.

In practice, a perfect reconstruction via a D/A converter is impossible to achieve as it relies on a perfect low-pass filter, described below, which cannot be constructed (Smith and Mersereau 1992).

Notationally, we have a continuous signal x_{t} that is sampled every T seconds to produce a sequence {x_{nT} }. For convenience, this will be referred to as x(n).

For any sequence x(n), the energy, E, is defined to be

E = ∑_{n=−∞}^{∞} x^{2}(n)

from Ludeman (1986).

Aliasing and the Nyquist rate

Sampling of continuous time signals has one major constraint associated with it, the problem of aliasing. This is where the reconstruction of a sampled signal is identical to the reconstruction of a signal with a different frequency. Figure 2.4 shows three signals with frequencies given by sin 2πt, sin 12πt and sin 8πt. The signals were all sampled with T = 0.2 and the reconstructions are shown. As can be seen, the sin 2πt signal’s reconstruction is very close to the original whilst the sin 12πt signal is aliased and appears to be a sin 2πt signal. The sin 8πt signal is also aliased but has the phase reversed.

Figure 2.5 indicates that the problem would be alleviated if the sample rate was higher, the question must be: what sample rate is high enough to prevent aliasing occurring? To justify the answer will involve using the Fourier transform which is detailed below — the mechanics are not important for this justification, merely the fact that the signal is moved into the frequency domain by the transform.

Figure 2.6 shows a continuous time signal with its Fourier transform and the sampled version with a matching Fourier transform. The triangles of the sampled signal’s Fourier transform show if the signal is aliased or not — if they overlap, aliasing has occurred. It is easy to see that there will be no aliasing if Ω_{0}T < π. This is equivalent to 1/T > Ω_{0}/π where Ω_{o} is the maximum finite radian frequency as shown in the diagram.

Now, Ω_{0} can be written as 2πf_{0} where f_{0} is a frequency in hertz, giving

\frac{1}{T} > 2f_{0}

Therefore, if the sampling rate 1/T is greater than 2f_{0}, there will be no

(a) x(t) = sin 2πt  \underrightarrow{sample}  _{(b) Unaliased}

(c) x(t) = sin 12πt  \underrightarrow{sample}  _{(d) Aliased}

(e) x(t) = sin 8πt

\underrightarrow{sample}

(f) Aliased with phase reversal

Figure 2.4: Various signals with their reconstructions.

aliasing of the signal and a perfect inversion to the original signal is possible. It follows that if a signal is sampled at a rate of n samples a second, the highest unaliased frequency in the reconstruction will be n/2 hertz.

The sampling rate 2f_{0} where f_{0} is the highest frequency present (or required) in a signal is known as the Nyquist rate.

Figure 2.5: Various signals and the sample positions when T = 0.2.

2.3 Signal Processing

As has been mentioned, signal processing is a wide and varied field with various mathematical transforms at the heart of most modern techniques. These transforms are realised in hardware, via filters, or mathematically, usually in software. As will be shown, there are comparatively few different types of mechanical filter, although there are many different transforms that make use of them.

Generally known transforms include the Fourier, Laplace, Walsh and z transforms although only the Fourier transform will be considered in any depth here. This is because it is one of a class of transforms that move from the time to the frequency domain and, as has been noted, this is highly desirable when sound data are being analysed. A more recent transform in this general class, the wavelet transform will also be introduced.

⟶
↓Sample  ↓

⟶

Figure 2.6: Fourier transform of a sequence resulting from sampling a continuous time signal.

2.3.1 Digital Filters

There are two main types of filters available today. The oldest, and most common, type are analogue filters that work on a continuous time signal whilst newer, digital filters work on discrete time signals.

The reasons for using digital filters, and hence the increasing usage of digital and sampled data, are many. Digital filters have been shown to be

“stable, predictable and repeatable. There is little degradation … and noise is limited …. The basic overall structure remains fixed yet allows changes via programming.” (Ludeman (1986), p17)

In addition, costly components can be replaced by software on a general purpose computer, decreasing marginal costs and design time. Analogue filters also suffer from degradation over time, even when not used, and short term drift through temperature and voltage fluctuations, causing further reliability problems.

Despite these advantages, digital filters have traditionally had two problems, the cost of the computer (or A/D, digital filter, D/A equipment chain) and the problem of sampling very high frequency signals. Recently however, these factors have been generally overcome by the major price/performance improvements seen in the computer industry, although the price/performance of the faster, broadband analogue filters is still some years ahead of their digital equivalents (Grant et al. 1989).

With cheap computing power and improved sampling hardware, digital filters are starting to replace analogue filters as the main method of signal processing, this has been further boosted by the appearance of dedicated DSP (Digital Signal Processing) chips being built in to computers of all types.

For all these reasons, the following discussion will concentrate on the digital filters to the exclusion of their analogue counterparts.

Types of Digital Filter

There are broadly two parameters for digital filters. The first indicates what frequency bands are passed through by the filter. The second indicates whether the filter is Infinite Impulse Response (IIR) or Finite Impulse Response (FIR). These two criteria will be looked at in turn.

Filters can be split into four main band types. These are low-pass, highpass, band-pass and band-stop filters. Figure 2.7 shows the ideal (solid line) and practical (dotted line) frequency responses of these filter types as frequency increases.

The ideal filters are not realisable mechanically, so actual filters have the transition stage shown by the dotted line where not all the energy in these frequencies are stopped. This transition ramp can be made steeper but can not be eliminated. Ideal filters are also known as brick-wall filters for obvious reasons.

Although there are four basic filter types, only the low-pass filter is ever

(a) Low-pass  (b) High-pass

(c) Band-pass  (d) Band-stop

Figure 2.7: Frequency responses of filters

actually needed as the other three types can all be built from low-pass filters. This has obviously meant that a lot of research has gone into achieving the best possible low-pass filter. As Ludeman (1986) shows, any low-pass, highpass, band-pass or band-stop filter can be obtained from transforms of a normalised low-pass filter — full details of these transforms are given in Table 3.2, page 128 of his text.

Digital filters all follow the basic form of producing a response sequence, y(n), from the inputs x(n) using the formula

y(n) = ∑_{k=0}^{M} c_{k}x(n − k) + ∑_{j=1}^{N} d_{j}y(n − j)

where the M + 1 coefficients c_{k} and the N coefficients d_{j} are fixed and designate the filter response with the filter producing every output value from the current and M previous input values and from its own N previous output values. Therefore, if N = 0, so that previous output values are not used, the filter is a finite impulse response filter and if N ≠ 0 the filter has an infinite impulse response (Press et al. 1992). Of course, an IIR filter does not necessarily have an infinitely long response, it is just capable of such.

The c_{k}’s and d_{j}’s are linked with the filter response function ϑ(f) which is given by

ϑ(f) = \frac{∑_{k=0}^{M} c_{k} exp (−2πik(fΔ))}{1 − ∑_{j=1}^{N} d_{j} exp (−2πij(fΔ))}

where Δ is the sampling interval. For FIR filters, the denominator is 1. The actual problem of filter design is based on the the inverse of this equation however. Given a desired filter response, the coefficients c_{k} and d_{j} need to be defined.

This area leads to a rich and varied literature with many texts being dedicated to solving this problem. This work is not about digital filter design, however fascinating this area might be, so instead, it is recommended that a suitable text, such as Ludeman (1986) is consulted.

An alternative to designing such a filter, although one which is not very appropriate for real-time applications, is to transform the data to the frequency domain before processing and then transform the altered data back to the time domain afterwards. The transforms to and from the frequency domain are performed by the Fourier transform which will now be discussed in greater detail.

2.3.2 The Fourier Transform

The Fourier transform is well known in mathematics and engineering, because it is highly useful to metamorphose a signal from the time to the frequency domain. This domain shifting allows certain features of signals to be examined and extracted more easily than would otherwise be the case.

To begin, the Fourier transform will be defined and some of its properties shown. Next, some methods to quickly compute the transform for special subsets of data will be given, and finally some applications of the transform will be introduced.

The basic principle of Fourier analysis is the representation of a set of data in terms of sinusoidal functions. One reason for choosing sinusoidal functions over and above other periodic functions is their simple behaviour under a change of time scale (Whalen 1971). A sinusoid with a period of 2π/ω may be expressed in the form

f(t) = R cos(ωt + ϕ)

where R is the amplitude and ϕ the phase. If the time variable is changed to u where u = (t − a)/b, say, to give a change of scale and of origin, we have

g(u) = f(a + bu)
= R cos(ωbu + ϕ + ωa)
= R′ cos(ω′u + ϕ′)

where R′ = R, ω′ = ωb and ϕ′ = ϕ + ωa, giving an unchanged amplitude, a frequency that is multiplied by the reciprocal of the change in the time scale, and a phase altered by the frequency of the sinusoid and the origin. As the time origin of a sinusoid is usually arbitrary, these relations are useful in time data analysis.

The Fourier transform will now be defined. Obviously, the transform could be defined by way of the Fourier series but it is simpler and more consistent to give the definition of the Fourier transform directly. Some notation is needed before proceeding:

Let ∥f∥_{1} be defined by

∥f∥_{1} = ∫_{−∞}^{∞} |f(x)|dx

Now, if the integral converges, then we say that ∥f∥_{1} is finite and takes the value to which the integral converges, otherwise (if the integral diverges), we say that ∥f∥_{1} = ∞.

The definition of the Fourier transform, f̂, can now be stated for a function f for which ∥f∥_{1} is finite (Walker 1991):

f̂(u) = ∫_{−∞}^{∞} f(x) exp(−i2πux)dx

The Fourier transform has four important properties; proofs are not given, but may be found in many of the standard texts on the Fourier transform (e.g. Nussbaumer 1981; Pettit 1985). The four properties are linearity, scaling, shifting and modulation.

Linearity implies that for all constants a and b, af +bg → af̂+bĝ. Scaling means that for all positive constants p, f(x/p) → pf̂(pu) and f(px) → (1/p)f̂(u/p). Shifting gives the result that for every real constant c, f(x − c) → f̂(u) exp(−i2πcu) and modulation means that f(x) exp(i2πcx) → f̂(u− c).

An important feature of the Fourier transform is that it can be inverted to move from the frequency domain back to the time domain. This allows a sound, for example, to be altered in the frequency domain and then reconstructed in the time domain and played back in order to hear the effect of the changes.

Again, proofs are available in many texts (e.g. Walker 1988) but the simplest inversion theorem is as follows (from Walker 1991):

If f is continuous and ∥f∥_{1} and ∥f̂∥_{1} are both finite, then

f(x) = ∫_{−∞}^{∞} f̂(u) exp(i2πux)du

for all x in ℜ^{2}. This integral is usually known as the inverse Fourier transform.

Finally, before looking at the discrete version of the transform, the total power^{1} of a signal is given by Parseval’s theorem:

Total Power = ∫_{−∞}^{∞} |f(t)|^{2}dt = ∫_{−∞}^{∞} |f̂(t)|^{2}dt

This means that the total power is the same irrespective of whether it is computed in the time or frequency domain (Press et al. 1992).

1power ∝ amplitude^{2}

The Discrete Fourier Transform

In this section, details of the Discrete Fourier Transform (DFT) will be given. The DFT is important as it can be used to approximate the Fourier transform, which is not in itself of great value for digital signal processing, and there are various methods, detailed below, to calculate it extremely quickly, which is highly useful.

The DFT is the equivalent of the Fourier transform for sampled digital signals — usually sampled continuous signals as is the case for the underwater sound data. In a nutshell, the DFT is an estimate of the Fourier transform from a finite number of its sampled points.

Suppose that we have N consecutive sampled points x_{t}, t = 0, 1, 2, …, N− 1 with a sample rate given by 1/T (as before). For ease of calculation, constrain N to be even. If the function x(t) is non-zero only in a finite time span, x_{T} should totally enclose this interval, otherwise x_{T} should be a good representation of what x(t) looks like at all times.

With N samples, we obviously cannot produce more than N transform values so estimates of the transform will be produced only at the discrete values f_{n} ≡ n/NT, n = −N/2, …, N/2 where the extreme values of n correspond to the lower and upper Nyquist limits. Although there are N + 1 numbers in this sequence, the extreme values are not independent (they are, in fact, equal) which gives the required N independent numbers.

We can now produce the Fourier transform by a discrete sum as shown by Press et al. (1992):

f(f_{n}) = ∫_{−∞}^{∞} f̂(u) exp(i2πuf_{n})du

= ∑_{t=0}^{N−1} x_{t} exp(i2πu_{t}f_{n})T

= T ∑_{t=0}^{N−1} x_{t} exp(\frac{i2πtn}{N})

This final summation is the discrete Fourier transform of the N points x_{t}. If we let H_{n} denote this, so that

H_{n} = ∑_{t=0}^{N−1} x_{t} exp (\frac{i2πtn}{N})

then the relationship between the DFT of a set of numbers and their continuous Fourier transform when they are viewed as a sample from a continuous function is

f(f_{n}) = T H_{n}

following Press et al. (1992).

It should be noted that although n is defined as varying from −N/2 to N/2, the same results hold with n in H_{n} in the range 0, …, N − 1 due to the fact that H_{n} is periodic in n with period N (which means that H_{−n} = H_{N−n}, n = 1, 2, …).

The formula for the inverse discrete Fourier transform is given by:

x_{t} = \frac{1}{N} ∑_{n=0}^{N−1} H_{n} exp (\frac{−i2πtn}{N})

which has the same form as the formula for H_{n} with a change of sign in the exponential, allowing a routine that calculates the DFT to also calculate the inverse with very few changes.

It should also be noted that the discrete form of Parseval’s theorem for the total power is

∑_{t=0}^{N−1} |x_{t}|^{2} = \frac{1}{N} ∑_{n=0}^{N−1} |H_{n}|^{2}

Methods

The discrete Fourier transform is, at first sight, an O(N^{2}) process with N^{2} complex multiplications and some additional computation. However, the DFT can be calculated using the Fast Fourier Transform (or FFT ) which reduces the process to one of O(N ln N) operations. This is a highly impressive reduction, for instance, for N = 1024 the FFT involves only 0.488% of the operations of the DFT which can reduce computation of such a DFT from weeks to seconds of computer time. The FFT became popularised in the mid-1960s with the work of Cooley and Tukey (1965). We now know (Ludeman 1986) that efficient methods for calculating the DFT have been discovered and implemented independently by several individuals since Gauss in 1805.

Whilst the FFT is a useful discovery, there is one potential drawback. Many of the algorithms for calculating it require the data length, N, to be an integer power of 2. At the simplest level, this means that data must be padded with zeros to expand the length to be a power of 2. There are, however, algorithms that allow any length of data to be Fourier transformed so this is not a practical problem.

Press et al. (1992) give a very clear derivation of the FFT, based on the work of Danielson and Lanczos in 1942 which will be reproduced here.

Danielson and Lanczos proved that a DFT of length N can be expressed as the sum of two Fourier transforms, each with length N/2. One is formed from the even points of the original values and the other from the odd points. The proof is:

F_{k} = ∑_{j=0}^{N−1} exp (\frac{2πijk}{N} ) f_{j}

= ∑_{j=0}^{N/2−1} exp (\frac{2πik(2j)}{N}) f_{2j} + ∑_{j=0}^{N/2−1} exp (\frac{2πik(2j + 1)}{N}) f_{2j+1}

= ∑_{j=0}^{N/2−1} exp (\frac{2πikj}{N/2} ) f_{2j} + W^{k} ∑_{j=0}^{N/2−1} exp (\frac{2πikj}{N/2} ) f_{2j+1}

= F_{k}^{e} + W^{k}F_{k}^{o}

where W ≡ exp(2πi/N), F_{k}^{e} denotes the kth component of the Fourier transform of length N/2 formed from the even components of the original f_{j}’s and F_{k}^{o} is the corresponding transform formed from the odd components.

This Danielson-Lanczos Lemma can be used recursively, so, having reduced the problem of computing F_{k} to that of computing F_{k}^{e} and F_{k}^{o}, we can do the same for F_{k}^{e} to produce F_{k}^{ee} and F_{k}^{eo}. As we assume that the original N is a power of 2 we can continue applying the lemma until all the data have been divided into transforms of length one. In fact, for every pattern of ln N e’s and o’s there exists a one-point transform that is just one of the input numbers f_{n}. So, for some n, F_{k}^{eoeeoeo...oee} = f_{n}

It is now necessary to identify which value of n corresponds to which pattern of e’s and o’s in the above equation. Luckily, this is easy; if the pattern of e’s and o’s is reversed and e set to 0 whilst o is set to 1 then the result is the binary value of n. It is this that makes computer processing of FFTs possible. Suppose we take the original vector of data f_{j} and rearrange it in bit-reversed order (so positions 0, 1, 2, 3 become 0, 2, 1, 3 for example) then the points as given are the one-point transforms. Combining adjacent pairs gives the two-point transforms, adjacent pairs of pairs give the fourpoint transforms, and so on until the two halves of the whole data set are combined into the final transform. Each combination takes O(N) operations and there are ln N combinations which gives an O(N ln N) operation. This algorithm is generally known as a decimation-in-time or Cooley-Tukey FFT. There is also a decimation-in-frequency or Sande-Tukey FFT algorithm which performs the ln N iterations on the input data first, and then rearranges the output values into bit-reverse order.

There are further refinements of the general FFT procedure which divide the initial data into sets of size 4 or 8 and optimised transforms performed on these small sets using special properties of that sized transform. These are known as base (or radix) 4 or 8 FFTs and can be 20–30% faster the general FFT shown above.

The FFT can be easily expanded to two or more dimensions which makes it suitable for image processing applications as well, indeed the JPEG (Joint Picture Expert Group) encoding standard is based on the fast Fourier transform to compress picture data. An overview of the JPEG algorithm and its suitability as a picture compression method is given in Lammi and Sarjakoski (1995).

The Power Spectrum

The power spectrum of a signal (also called the power spectral density — PSD) is a measure of how the total power of a signal is spread between the various frequencies it contains. There are many ways to measure this for a sampled signal so the periodogram will be looked at as this, historically, is the most widely used estimate of the power spectrum. If we have an N-point sample of the function x(t) and we compute its discrete Fourier transform

X_{k} = ∑_{j=0}^{N−1} x_{j} exp (\frac{2πijk}{N} ) k = 0, …, N − 1

then the periodogram is defined at \frac{N}{2} + 1 frequencies to be

P (0) = P (x_{0}) = \frac{1}{N^{2}} |X_{0}|^{2}

P (x_{k}) = \frac{1}{N^{2}} (|X_{k}|^{2} + |X_{N−k}|^{2}) k = 1, 2, …, (N/2 − 1)

P (x_{c}) = P (x_{N/2}) = \frac{1}{N^{2}} |X_{N/2}|^{2}

where x_{c} is the Nyquist frequency and x_{k} is defined only for the zero and positive frequencies

x_{k} ≡ \frac{k}{NT} = 2x_{c} \frac{k}{N} k = 0, 1, …, N/2

This, then, is the periodogram. However, it is obviously not a perfect estimator of the true power spectrum of the signal, if only because frequencies are “binned” so a P (x_{k}) value will be some kind of average power for a range of frequencies. The expectation is that the P (x_{k})’s are some kind of average of P (x) over a narrow window function centred on the x_{k}’s. For the periodogram, this window function, given as a function of s the frequency offset in bins, is

W (s) = \frac{1}{N^{2}} (\frac{sin(πs)}{sin(πs/N)})^{2}

It should be noted that W (s) falls off at a rate of approximately (πs)^{2} which is not very rapid. This results in leakage between frequencies in the periodogram estimate. The solution to this problem is data windowing which will now be discussed.

Data windowing has one goal, to minimise the leakage in bin k from frequency components that are s bins away. Now, with no data windowing, when we process a run of N sampled points we are actually multiplying an infinite run of sampled data x_{t} by a window function in time that is zero except during the total sampling time NT when it is unity. This is a square window and the reason it has large amounts of leakage is that it switched on and off instantaneously. To overcome the leakage effect, a window that changes gradually from zero to a maximum and back again should be used, with the proviso that it must still be zero outside the sampling time NT of course.

With such a window, w_{t}, the periodogram equations are given by:

P (0) = P (x_{0}) = \frac{1}{W_{ss}} |D_{0}|^{2}

P (x_{k}) = \frac{1}{W_{ss}} (|D_{k}|^{2} + |D_{N−k}|^{2}) k = 1, 2, …, (\frac{N}{2} − 1)

P (x_{c}) = P (x_{N/2}) = \frac{1}{W_{ss}} |D_{N/2}|^{2}

where D_{k} ≡ ∑_{t=0}^{N−1} x_{t}w_{t} exp(2πitk/N), k = 0, …, N − 1 and W_{ss} is window squared and summed, W_{ss} ≡ N ∑_{t=0}^{N−1} w_{t}^{2}. Obviously, x_{t} is as given above. This gives the general form of W (s) to be

\begin{aligned}W (s) & = \frac{1}{W_{ss}} | ∑_{k=0}^{N−1} exp(2πisk/N)w_{k}|^{2} \\ & ≈ \frac{1}{W_{ss}} | ∫_{−N/2}^{N/2} cos(2πsk/N)w(k − N/2)dk|^{2} \end{aligned}

The approximation is useful for practical estimates when s << N and the window is left-right symmetric (as most are). The function w(k − N/2) is a smooth, continuous function passing through the points w_{t} (Press et al. 1992).

Some of the more common windows functions are the rectangular, Welch, Bartlett and Hanning windows. For most applications there is no practical differences between them, the Welch window is usually a good compromise choice. Figure 2.8 shows these typical windows.

(a) Rectangular  (b) Hanning

(c) Welch  (d) Bartlett

Figure 2.8: Window functions commonly used in FFT power spectrum estimation.

The Welch window is given by

w_{t} = 1 − (\frac{t − \frac{1}{2}N}{\frac{1}{2}N} )^{2}

So far, this discussion has concentrated on looking at a single periodogram for the sample data. Many applications have signals that change over time however, so an adaption of the periodogram was developed to have a time adaptive power spectrum. This technique is know as the spectrogram. The procedure is simple, for a data set of M points, sets of N points (where N is a power of 2) are taken and a periodogram produced for the sets of N points using a suitable window. The N point set moves along the data set until all the M points have had periodograms produced.

The N point set can be moved along either in blocks of N points so that each point is used for one periodogram or in blocks of less than N points so that blocks overlap.

The results are displayed as in Figure 2.9. Each column of pixels is one periodogram with frequency on the vertical axis and time on the horizontal. The values of the periodogram determine the colour of the pixels, with dark pixels showing low power levels and light pixels being high power.

Figure 2.9: An example spectrogram — the word “operator”.

As can be seen, the high power areas are in a relative narrow band, although there are harmonics at higher frequencies. These spectrograms will be used in later chapters to give visual impressions of the data used.

Although the Fourier transform and the derived periodograms are the most common method of viewing and analysing sound data, the technique does have some problems which are hard to overcome. These problems, and a different transform to overcome them, are detailed in the next section, the wavelet transform.

Signal to Noise Ratio

A brief note should be made here about the signal to noise ratio (SNR). This is a standard measure of how “clean” the signal is, i.e. how much noise it contains.

The SNR is measured in decibels (dB) and is defined by:

SNR = 10 log_{10} (\frac{S}{N} )

where S is the power of the signal (where S will be proportional to the squared amplitude of the signal component of the sound) and N is the power of the noise (proportional to the squared amplitude of the noise component). However, the SNR is only valid for statistically stationary noise. In the case of underwater sound this is not the case (as will be shown in Chapter 3) but this measure is very widely used and therefore will be used throughout this work, following the precedent set by Smith (1993).

Now, as the power of a sound is proportional to the square of the amplitude of the sound then, for a sound which is all noise and no signal, the Noise Level (NL), an similar measure to the SNR, is given by:

NL = 10 log_{10} (\frac{1}{σ^{2}} )

where σ^{2} is the variance of the power of the noise.

2.3.3 The Wavelet Transform

Fourier transforms are the ideal decompositions for stationary signals. Stationary signals are ones whose properties are statistically invariant over time, such as white noise. Any unexpected events that appear must have their probabilities known in advance (Meyer 1993b). Our signals, as discussed in the introductory chapter, have transients appearing at unknown times for unknown durations and a background that changes over time with no readily discernible sequence, which makes our signals non-stationary. It is for this reason that we look for techniques other than Fourier analysis, as problems occur when trying to adapt Fourier analysis to signals it was not designed to decompose with any degree of success.

In signal processing, it is advantageous to look at time-scale decompositions as opposed to the Fourier frequency decomposition. In some situations (where, for example, one has a stationary background with unpredictable transients) it may be advantageous to look at time-frequency decompositions, which can be performed by a special type of wavelet know as Malvar wavelets, but in our situation, where even the background is non-stationary, the most suitable technique is time-scale decomposition.

As the name implies, time-scale analysis involves looking at a signal in time and scale, where scale is analogous to the concept of scale in cartography. At a given scale, the signal is replaced by the most detailed approximation possible, so that, as the scale increases, the approximation becomes more exact.

This obviously requires a range of scales to be looked at, with the analysis involving an examination of the changes as the scale changes. This technique of examining increasingly finer scale levels is known as multi-resolution analysis and is equivalent to an atomic decomposition where the atoms are Grossmann-Morlet wavelets (Meyer 1993b).

These wavelets are defined around the concept of a mother wavelet, ψ(t) of the real variable t. It has to be well localised and oscillating (resembling a wave). The localisation is that ψ(t) decreases rapidly to zero as |t| tends to infinity. Oscillation suggests that the mother wavelet vibrates like a wave, so we require that the integral of ψ(t) be zero and that the same holds true for the first m moments of ψ (Meyer 1993b). This can be expressed as

0 = ∫_{−∞}^{∞} ψ(t)dt = ⋯ = ∫_{−∞}^{∞} t^{m−1}ψ(t)dt

For a function ψ(t) to be a mother wavelet, there are other conditions to be satisfied, mostly to allow that wavelet family to be used on function spaces other than L^{2}(ℜ), but these will be detailed below. For now we will assume that we have a mother wavelet that fulfils all the necessary conditions and the construction of the wavelet family will be demonstrated.

The mother wavelet is used to construct the other wavelets in the family by means of two processes, dilation and translation. Dilation alters the width of the wavelet to give the different scales and translation moves the wavelet across the signal. So, if ψ(t) is the mother wavelet, then the wavelets, ψ_{j,k}(t), are given by

ψ_{j,k}(t) = 2^{\frac{j}{2}} ψ(2^{jt} − k) j > 0, 0 ≤ k < 2^{j}

where j represents the dilation number and k is the translation number. 2^{j/2} is a scale factor to normalise ψ_{j,k} so that ∥ψ_{j,k}∥ = ∥ψ∥.

Now, if ψ is a mother wavelet, then the set of functions ψ_{j,k} form an orthonormal basis and it is therefore possible to use the wavelets as the basis for an expansion of a function,f(·), of interest. The wavelet series representation of such a function is then

f(t) = ∑_{j,k} d_{j,k}ψ_{j,k}(t),

where the wavelet coefficients, d_{j,k}, are given by

d_{j,k} = ∫_{−∞}^{∞} f(t)ψ_{j,k}(t)dt =< f, ψ_{j,k} >

where < ., . > denotes the inner product (Nason and Silverman 1994).

Although Grossmann and Morlet showed this and then stated that this new technique was suitable for the analysis of transient signals, they were derided by established signal processing experts. However, Stephane Mallat showed that this technique had been in use for several years under different names such as quadrature mirror filters and pyramid algorithms (Mallat 1989). He created a method whereby this signal analysis technique was shown to be in use with differing presentations and this synthesis spurred on many new developments, including Daubechies’ discovery of a new system of mother wavelets which will be detailed later.

The oldest example of a function ψ for which the ψ_{j,k} constitute an orthonormal basis for a function in L^{2}(ℜ) is the Haar function defined by ψ(t) = 1 on [0, 0.5), ψ(t) = −1 on [0.5, 1] and ψ(t) = 0 elsewhere. This function at the first two levels of detail is shown in figure 2.10 (after Strang 1993).

Figure 2.10: Haar wavelet (ψ(t)) and the next detail level.

The Haar wavelet was developed in 1909 (although it was not called a wavelet until much more recently). It is now recognised as a special case amongst Daubechies’ wavelets.

As stated above, for a function to be classed as a mother wavelet, extra conditions need to be imposed above and beyond the localisation and oscillation conditions described above. Meyer (1993b) describes these extra conditions as

the regularity of the wavelet ψ

the decay at infinity of ψ, and

the number of vanishing moments.

Depending which of these conditions is fulfilled and exactly how they do so, produces the various types of wavelets available in the general time-scale decomposition field. These various types of wavelets will be discussed next.

Types of Wavelet

Although the earliest wavelet system was, as we have seen, developed by Haar in 1909 (a mere 100 years after Fourier produced the Fourier system), there are two main criticisms which took another 80 years to overcome. The first criticism is that although the Haar system attempts to model a continuous function, the wavelets are not themselves continuous which gives a lack coherence. The other, more serious, problem is that only continuous functions which are square-integrable on [0, 1] can be modelled successfully (Meyer 1993b).

To overcome these objections, the three properties listed above need to be added to the definition of a mother wavelet. Following on from the work of Mallat, Daubechies was able to complete Haar’s original work by requiring that ψ(t) has m continuous derivatives, compact support and that the first r + 1 moments vanish. She was thus able to construct, for each integer r, an orthonormal basis for L^{2}(ℜ) of the form

2^{\frac{j}{2}} ψ_{r}(2^{jt} − k)

which has support of ψ_{r} on [0, 2r + 1], ψ_{r}(t) has γr continuous derivatives (where γ is about one fifth) and ∫_{−∞}^{∞} t^{r}ψ_{r}(t)dt = 0 ∀r. When r = 0 this gives the Haar basis (Daubechies 1992).

One of the major achievements of the wavelet basis is that both orthogonality and local support are achieved. Orthogonality implies that

∫ ψ(x)ψ(2^{jx} − k)dx = 0

for all integers j and k where j is dilation and k is translation. Thus ∑ b_{jk}ψ_{jk}(x) is a multi-resolution of f(x) in which b_{jk} carries information pertaining to f near frequency 2^{j} at x = 2^{−j}k (Strang 1989). Such orthogonality with local support is not easy to achieve, a sine wave truncated at 0 and 2π, for instance, is not orthogonal to all its dilations. In fact, this dual constraint leads to irregular basis functions — all calculations are performed recursively.

The scaling function, ϕ is constructed through a dilation equation

ϕ(x) = ∑ c_{k}ϕ(2x − k)

where a solution normalised by ∫ ϕdx = 1 is sought. The initial requirement on the coefficients c_{k} comes from multiplying by two and integrating to give

2 ∫ ϕdx = ∑ c_{k} ∫ ϕ(2x − k)d(2x − k) yields ∑ c_{k} = 2.

∑ c_{k} = 2 ensures uniqueness of ϕ, however a smooth solution is not ensured. A method of constructing the ϕ will now be given (after Strang (1989), construction 3). A wavelet ψ(x) is a further combination of the translates ϕ(2x − k) involving the same c_{k}.

Suppose that ϕ is known at the integers x = j. The basic dilation equation recursion gives ϕ at the half-integers then at the quarter-integers and ultimately at all points x = k/2^{j}. This is very quick and is the method used in most practical wavelet calculation routines.

The values of ϕ at the integers come from an eigenvector. For example, for Daubechies’ extremal phase wavelet of order 4 (shown in figure 2.11) the c_{k}’s are 1/4(1+\sqrt{3}), 1/4(3+\sqrt{3}), 1/4(3−\sqrt{3}) and 1/4(1−\sqrt{3}). Using these and setting x = 1 and x = 2 in the dilation equation whilst using the fact that ϕ = 0 unless 0 < x < 3 we have:

\begin{aligned}ϕ(1) & = 1/4(3 + \sqrt{3})ϕ(1) + 1/4(1 + \sqrt{3})ϕ(2) \\ ϕ(2) & = 1/4(1 − \sqrt{3})ϕ(1) + 1/4(3 − \sqrt{3})ϕ(2) \end{aligned}

This is ϕ = Lϕ with L_{ij} = c_{2i−j}. The eigenvalues are 1 and 1/2 and the eigenvector for λ = 1 has components ϕ(1) = 1/2(1 + \sqrt{3}) and ϕ(2) = 1/2(1 − \sqrt{3}) which are the extreme heights of figure 2.11.

In fact, from ϕ(1) and ϕ(2), recursion produces all that is needed.

Figure 2.11: Daubechies’ Extremal Phase order 4 wavelet.

It would be pointless to reproduce all the mathematics behind the following properties that are a result of this approach to finding ϕ. Instead, the interested reader is advised to pursue some or all of the following texts: Strang (1989), Daubechies (1993), Daubechies (1990), Meyer (1993b) and, of course, Daubechies (1992).

The important point is that scaling functions produced from such dilation equations produce wavelets with the following properties:

• Smooth functions can be approximated with error O(h^{p}) at every scale h = 2^{−j}

• The first p moments of the wavelet ψ(x) are zero, i.e.:

∫ x^{m}ψ(x)dx = 0 for m = 0, …, p − 1

• The recursion matrix L_{N} that determines ϕ at the integers has the eigenvalues 1, 1/2, …, (1/2)^{p−1}.

This last property is very important as it is this which makes ϕ(x) smoother as p increases so that higher order wavelets are much smoother than the low order wavelet shown in Figure 2.11.

We can now define a wavelet, ψ(x) using the scaling function and the c_{k}’s:

ψ(x) = ∑ (−1)^{k}c_{1−k}ϕ(2x − k).

Daubechies great feat was to construct a family of such wavelets which incorporated Haar’s wavelet as the simplest case. The wavelets have accuracies 1, 2, 3, … having 2, 4, 6, … coefficients c_{k} not equal to zero. As the accuracy increases so does the smoothness. Compare Figure 2.12, an extremal phase order 10 wavelet, with the same families order 4 wavelet shown in Figure 2.11.

Figure 2.12: Daubechies’ Extremal Phase order 10 wavelet.

It should also be noted that the wavelet transform can be extended to higher dimensions with ease to allow the processing of images and other multi-dimensional data. Daubechies (1992) gives details but as sound is being analysed in this work, the one-dimensional form of the transform will be the only one considered here.

The Discrete Wavelet Transform

Just as with the Fourier transform, there is a discrete form of the wavelet transform, the Discrete Wavelet Transform (DWT). As well as allowing digital data to be subjected to a wavelet transform, this also has a fast version know as the Fast Wavelet Transform or FWT much like the FFT. The wavelet form has the advantage of involving even less calculation than the Fourier equivalent.

The DWT is a pyramid algorithm as developed by Mallat (1989). The implementation used here was developed by Nason and Silverman (1994) for S-Plus and the following description of the wavelet transform algorithm is based on the description given in Section 7 of that paper.

It is assumed that a scaling function, ϕ(x), has been obtained by the method described above and there is thus a mother wavelet, ψ(x), which has been used to obtain a wavelet family by the equation

ψ_{jk}(x) = 2^{\frac{j}{2}} ψ(2^{j}x − k)

and that the wavelet coefficients are denoted by the vector d_{k}^{j} with j being the scaling and k the translation.

Following Mallat (1989), the process starts with a set of data c_{0}^{M}, …, c_{N−1}^{M} where the superscript is the resolution level and the subscript is the coefficient within that level. The superscript M represents data at the highest resolution level and there are N items of data. The transform has M = ln N stages with each such stage producing a set of smoothed coefficients c and detail coefficients d. Figure 2.13 shows some of the steps for a DWT on 16 items of data. Hopefully this will be a useful guide during the description that follows.

The algorithm involves applying the wavelet coefficient matrix described above in a hierarchical fashion, first to the full set of data points and then to the “smooth” vector which is half the length of the previous vector. The final outcome contains the remaining “smooth” component and the “detail” components from each level. Increasing the initial number of data points to a higher power of two would simply increase the number of transformation

[\begin{matrix}x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \\ x_{5} \\ x_{6} \\ x_{7} \\ x_{8} \\ x_{9} \\ x_{10} \\ x_{11} \\ x_{12} \\ x_{13} \\ x_{14} \\ x_{15} \\ x_{16} \end{matrix}] ⟶^{transform} [\begin{matrix}s_{1} \\ d_{1} \\ s_{2} \\ d_{2} \\ s_{3} \\ d_{3} \\ s_{4} \\ d_{4} \\ s_{5} \\ d_{5} \\ s_{6} \\ d_{6} \\ s_{7} \\ d_{7} \\ s_{8} \\ d_{8} \end{matrix}] ⟶^{permute} [\begin{matrix}s_{1} \\ s_{2} \\ s_{3} \\ s_{4} \\ s_{5} \\ s_{6} \\ s_{7} \\ s_{8} \\ d_{1} \\ d_{2} \\ d_{3} \\ d_{4} \\ d_{5} \\ d_{6} \\ d_{7} \\ d_{8} \end{matrix}] ⟶^{transform} [\begin{matrix}S_{1} \\ D_{1} \\ S_{2} \\ D_{2} \\ S_{3} \\ D_{3} \\ S_{4} \\ D_{4} \\ d_{1} \\ d_{2} \\ d_{3} \\ d_{4} \\ d_{5} \\ d_{6} \\ d_{7} \\ d_{8} \end{matrix}] ⟶^{permute} [\begin{matrix}S_{1} \\ S_{2} \\ S_{3} \\ S_{4} \\ D_{1} \\ D_{2} \\ D_{3} \\ D_{4} \\ d_{1} \\ d_{2} \\ d_{3} \\ d_{4} \\ d_{5} \\ d_{6} \\ d_{7} \\ d_{8} \end{matrix}] ⟶^{etc} [\begin{matrix}S_{1} \\ D_{1} \\ D_{1} \\ D_{2} \\ D_{1} \\ D_{2} \\ D_{3} \\ D_{4} \\ d_{1} \\ d_{2} \\ d_{3} \\ d_{4} \\ d_{5} \\ d_{6} \\ d_{7} \\ d_{8} \end{matrix}]

Figure 2.13: The Discrete Wavelet Transform for 16 data points.

and permutation stages. The final vector however, would always be a hierarchy of detail coefficients and the smooth coefficients. As each stage is an orthogonal linear operation, the whole process is also an orthogonal linear operator (Press et al. 1992).

It should also be noted that the DWT is, in effect, a series of iterative applications of two FIR digital filters. Firstly a “smoothing” (low-pass) filter to retrieve the “smooth” coefficients, c, and then a band-pass filter to retrieve the “detail” coefficients, d. These filters, h and g say, are known as quadrature mirror filters and are related by g(n) = (−1)^{n}h(1 − n) as shown by Daubechies (1990). It is thus only necessary to know one of the pair for any wavelet family. Also, the combined lengths of the new c’s and d’s is identical to the length of the old c’s.

One such operation will now be looked at in detail, decomposing level j into level j − 1. This decomposition will extract both a smoother signal at a lower resolution and the signal detail at that resolution.

The smoothing of the data at level j to the next level j − 1, is achieved via the smoothing filter h(n), 0 ≤ n ≤ N_{h} −1 with N_{k} nonzero coefficients. This smoothing is a convolution preceding a dyadic decimation (Mallat 1989):

c_{k}^{j−1} = ∑_{n∈Z} h(n − 2k)c_{n}^{j} = ∑_{m∈Z} h(m)c_{m+2k}^{j} = ∑_{m=0}^{m=N_{k}−1} h(m)c_{m+2k}^{j}

Nason and Silverman (1994) show the various inequalities which allow the coefficient indices, k to be computed but the above equation gives the smoothed data. To find the detail, the formula

d_{k}^{j−1} = ∑_{n∈Z} g(n − 2k)c_{n}^{j}

is used. Combining this with the equation above relating the filters g and h, we have

d_{k}^{j−1} = ∑_{n∈Z} (−1)^{n}h(2k + 1 − n)c_{n}^{j} = ∑_{m=0}^{m=N_{k}−1} (−1)^{1−m}h(m)c_{2k+1−m}^{j}

The end result of this algorithm is, as shown in Figure 2.13, a vector

DWT\{c^{M}\} = \{c_{0}^{0}, d_{0}^{0}, d_{0}^{1}, d_{1}^{1}, …, d_{0}^{M−1}, …, d_{\frac{N}{2} −1}^{M−1}\}

As can be seen, all coefficients except the first are detail values whilst the single smooth value, c_{0}^{0}, represents some mean of the data. For the Haar wavelets, this value is actually N times the sample mean.

The inverse DWT (or IDWT ) begins with the lowest resolution and builds up to full resolution. Using c_{0}^{0} and d_{0}^{0} the algorithm produces c_{0}^{1} and c_{1}^{1}. These c and d_{0}^{1} with d_{1}^{1} would then produce the next level until the final c_{k}^{M−1} are produced.

The general step, from level j to j + 1 is given by:

c_{n}^{j+1} = ∑_{k∈Z} h(n − 2k)c_{k}^{j} + ∑_{k∈Z} g(n − 2k)d_{k}^{j}

where g and h are the filters given above. Writing g in terms of h and limiting the summations to the possible values gives

c_{n}^{j+1} = ∑_{2k≥n+1−N_{k}}^{2k≤n} h(n − 2k)c_{k}^{j} + ∑_{2k≥n−1}^{2k≤N_{k}+n−2} h(1 + 2k − n)d_{k}^{j}

This needs no other knowledge than that used for the decomposition, so inverting it is just as fast.

Similar descriptions of this process are also given by, amongst others, Strang (1989), Daubechies (1990), Bock (1992) and Bock and Pliego (1992).

The only other point is how the boundary points are dealt with, usually the choices are either to treat the segment being decomposed as periodic or to make the end points symmetrical outside the segment of interest. Nason and Silverman (1994) give details on the difference these two approaches make to the wavelet transform.

The Fast Wavelet Transform, which involves a factorisation into very sparse matrices, takes the operation count to O(n). It involves log_{2} n steps — just like the FFT — but whereas each step of the FFT requires n/2 multiplications, the cost of each successive FWT step is cut in half by using the pyramid scheme created by Mallet. The total cost has a factor 1 + 1/2 + 1/4 + ⋯ that stays below 2 which is why the FWT is O(n). This makes it faster than the FFT which is O(n log n).

Comparison with Fourier Transform

At first sight, Fourier and wavelet transforms are very similar. Both decompose a signal based on basis functions (cos(·) and sin(·) in Fourier analysis, mother wavelets in wavelet analysis) and both move from the time domain to a frequency or frequency-like domain. This being the case, why should one be chosen over the other for any particular problem, and indeed, is there any evidence that one or the other is better for a class of signals?

The first contest between Fourier and wavelet methods was the MPEG (Motion Picture Experts Group) competition in Japan during 1991. This was designed to find the best way of compressing moving images into a certain size byte stream (150 kB/s), where “best” meant the highest visual quality according to the expert panel. Although the winner was a Fourier based technique, wavelets were still in their infancy so this was not too surprising (Strang 1993).

A later competition between the two was the Federal Bureau of Investigations (FBI) fingerprint database. The FBI has a fingerprint database in excess of 30 million records and the need for quick matches is extremely high. The obvious solution was to digitise all the records. To perform matches, the details of the ridge endings and bifurcations need to be retained (from a noisy background) and at 500 pixels per inch with 256 levels of gray scale, one record has over 10^{7} bytes of data. Compression of at least 20:1 was desired (Hopper and Preston 1992). When using the JPEG compression system, which is Fourier based, the ridges were broken and the image was unusable. With a wavelet based system however, the ridges were retained and the required compression ratios achieved. Such a system is now in use by the FBI for storage purposes (Meyer 1993a). However, recent research by Kidd (1995) has shown that with a suitably optimised quantisation matrix, JPEG results are at least on a par with wavelet based methods. As he admits however, if the wavelet based methods were also optimised, they could regain the lead. This is doubtless a rich area for future research.

The fingerprint data is certainly non-stationary and it is interesting to observe that Fourier analysis is not the best performer, but why should this be? There are several reasons but the most important ones are windowing and speed.

The speed improvement is not the main reason why wavelet transforms perform better — indeed, until the advent of hardware wavelet chips the FFT will be faster in practice — but it is a useful bonus. The FWT is a factor of log N faster than the FFT as described above. This is a useful increase on the already quick FFT.

The main advantage of the FWT however, is wrapped up in the necessity of windowing the data for the transforms. From the definition of the Fourier transform, we can see that if we require to discover information about any particular frequency, we must integrate over the whole time domain from t = −∞ to t = ∞. In many cases, we only require to know about the frequency behaviour over a short period of time. Indeed, in a real-time situation, with new data arriving constantly, there is no choice in the matter.

With the Fourier transform therefore, data must be windowed to produce the necessary time periods — various windows are described in the section entitled “The Power Spectrum” above. A windowed Fourier transform is known as the Short Time Fourier Transform (STFT). Such a transform does contain certain problems, as explained in Powell et al. (1995).

Firstly, the transform inherently loses time information since in integrates over time. To a large extent, the optimum size for the window will be dependent upon the data which could well be unknown in advance. If the time window employed is too large, a STFT may provide a good broad description, but lose individual events of short duration, such events occurring frequently in underwater sound data. Equally, if the time is too short then the overall picture is lost and long term events may be discarded. Window size, therefore, becomes an important factor in the analysis and results may be heavily dependent on the choice of this parameter. In general, it will not be possible to determine one window size which performs well across the entire range of sounds of interests, so differing sizes are needed in different sections. Moreover, for each different window size a complete STFT computation is required, since one cannot construct the STFT for smaller windows directly from one using a larger window or vice versa.

A second problem with using Fourier analysis, particularly relevant in the discrimination context, is that power spectra of recordings made in different circumstances must be levelled before the Fourier decompositions can be compared. Without this levelling of spectra, results are biased towards certain frequencies, and since the background characteristics of different types of underwater sounds are quite different — dependent on where and when they are recorded — serious comparison of the decompositions of different sounds cannot sensibly be undertaken without this additional, time consuming step.

Wavelets, however, provide time-frequency localisation, that is they can obtain information about frequency local to a given time point. From the description above, we can see that the coefficient f_{jk} gives information about the signal around time point 2^{−j}k at a frequency proportional to 2^{j}. This permits the wavelet transform to better “zoom in” on the irregularities of the signal and characterise them locally. It is also invariant to the length of the successive time windows which are decomposed; one can immediately obtain the wavelet decomposition of any required sub-window of time from that of a larger window without the need for further computation. Furthermore, there is no need to employ additional, usually ad-hoc, methods of levelling the power spectra; this is effectively performed automatically within the wavelet decomposition.

Thus, although the DWT is, in many ways, similar to the STFT, it overcomes the two main objections described, is faster and is still better able to describe short-lived, high frequency phenomena, such as transients in signals.

2.4 Reasons for Using Wavelet Analysis

The work in the remainder of this thesis is based on the use of wavelet analysis rather than Fourier analysis. There are several reasons for this approach and some will now be detailed.

The first is that when this work was started there was very little evidence of practical results obtained using wavelets. Most research tended to be on theoretical properties rather than an application based approach, this was doubly true when unusual data sets are concerned. To a great extent this is still the case so there is a need for some practical work to get a feel for how useful wavelets actually are in different situations. Fourier techniques have been exhaustively investigated however, making them a far less fruitful area for further research.

Secondly, as has been stated above, wavelet based techniques should be as good, if not better, than Fourier techniques with underwater sound data. It is an interesting exercise to see if the theoretical advantages of wavelets materialise in practice when real data are used.

There is one area where Fourier techniques do have a noticeable advantage and that is in the amount of software available that uses the FFT for analysis. Wavelet transforms are only just beginning to be explored which does limit the amount of software available unless one is willing to develop such software oneself — as was done by Nason (1993). Luckily there is enough software available to make wavelet analysis possible. Most of the software does limit one to using the Daubechies wavelets however, hence these are the wavelets used for the analysis presented here. The software used in this work is a combination of WaveThresh by Nason (1993) and S+Wavelets by Bruce and Gao (1994), both of which are implemented as plug-ins to the S-Plus program. Similar additions are also available for other statistical software such as Matlab and stand-alone wavelet packages are also now being produced.

Finally, concentrating on wavelet analysis has allowed the more common wavelet techniques to be “fine-tuned” towards our particular class of data, which is both interesting and innovative.

For all these reasons, the vast majority of the work done in this thesis is concerned with wavelet approaches. It is not a primary objective of the thesis to compare such approaches with those based on alternative transforms, but rather to investigate how wavelets can be moulded to our particular application and what problems are involved in that approach. Given these objectives, no emphasis is given in subsequent chapters to comparisons of how well Fourier techniques would or would not perform in all the particular areas investigated.

About the Data

3.1 Introduction

Having discussed the form of analysis that will be used, a fuller description of the data will now be provided. As has already been noted, the data were digitally sampled from an original analogue recording.

There are several different types of sound used in subsequent chapters; four real-life underwater sounds — dolphins, shrimps, seals and rain — and an artificial data set used for calibrating the procedures and for conducting comparative tests.

This chapter will consist of three main sections: a description of the recording circumstances and the data format, along with some statistical description of the theoretical distribution of this data; a second section detailing the features that can be expected in the true underwater sounds that need to be discovered, including some descriptive analysis of the four types of real data; and finally, a descriptive analysis of the artificial sound with some notes about its usage.

3.2 Technical Aspects of Underwater Sound

Data

This section concentrates on the technical aspects of underwater sound recording and analysis. After a description of the actual data format, some statistical theory will be presented to describe the underwater noise.

3.2.1 Recording Circumstances and Data Format

The four sets of real data, dolphins, shrimps, seals and rain, were all recorded in the same manner, although at differing times. A hydrophone, or underwater microphone, was suspended several meters under the surface of the sea (sometimes by a cord, at others by a metal rod) to record the aquatic sounds. This was attached to a high quality music cassette machine, similar to those found in domestic hi-fi systems. The cassette machine has a recording band of 0–16 kHz which is similar to the limits of human hearing Smith (1993).

In order to analyse these recordings on a computer, they must first be transferred from analogue to digital recordings by way of D/A converter, as described in Chapter 2. The Nyquist theory presented therein shows that in order to preserve all the recorded sound, a sampling rate of at least double the highest rate is needed. In fact, in order to avoid errors, an even higher rate is to be preferred, whilst still limiting subsequent analysis to the 0–16 kHz band.

For the data presented here, a sampling rate of 40.96 kHz, 2.56 times the highest frequency, has been adopted — this is the rate suggested by Kay Elemetrics Corporation in their Sona-Graph spectrum analyser which is used by the Defence Research Agency for such work. The actual digital storage of each sample is accomplished in 16 bits giving a range of -32768–32767 for each data point. This system is comparable to compact disc which has a sampling rate of 44 kHz and 16 bit recording.

The fidelity of this system easily retains all the features of the original analogue recording and is, in fact, beyond that of the original cassette used in the recording process.

3.2.2 Statistical Description of Underwater Sound

Several texts, such as Creasey et al. (1989a), Creasey et al. (1989b) and Urick (1975), have provided theoretical results for the probability densities of underwater sounds, whilst Smith (1993) gives an easy to follow summary.

The important points to arise from these papers are that underwater noise can be assumed to be formed from the combination of a large number of independent sources, meaning that the resulting noise will follow a normal distribution (by the central limit theorem) with a magnitude that has a Rayleigh PDF (Probability Density Function) and a power that has a negative exponential PDF. In cases where a sound is present, the magnitude follows a Rician distribution. The full derivation of these results can be obtained from Smith (1993).

As a point of interest, the probability density functions of the Rician and Rayleigh distributions are given by:

P_{Rician}(η) = \frac{η exp −(η^{2} + γ^{2})/2σ^{2}}{σ^{2}} I_{o}(\frac{γη}{σ^{2}} )

and

P_{Rayleigh}(η) = \frac{η exp −(η^{2})/2σ^{2}}{σ^{2}}

where, for underwater sounds, η is the magnitude of the sound and γ is the peak value of any signal present.

Despite these results, the problem of identification is not as easy as it first appears. It seems as if one merely needs to determine whether the magnitude of the incoming sound has changed from a Rayleigh to a Rician distribution; if so then a signal has been determined. The next section demonstrates that this is not such an easy problem, as some signals should, for practical purposes, be treated as noise, thus invalidating this method of analysis.

3.3 The Features of Underwater Sound

In this section, the various types of underwater sounds are described in some detail; these are recordings of dolphins, shrimps, seals and rain. The latter is used as a “control”, for determining whether analytical methods are behaving in a suitable manner — with rain, no signals should ever be detected. To begin though, some of the inherent features and problems with underwater sound analysis are described.

3.3.1 Features Inherent in Underwater Sound

The hydrophone used to record these sounds responds to changes in pressure, as that is the basis of sound transmission. We wish to split the incoming sound into signal or noise but given that the hydrophone picks up all pressure changes this presents a problem.

For example, suppose that an underwater volcano explodes some distance from the hydrophone. This will produce a large change in pressure for some time and will have the effect of increasing the magnitude of the incoming sound tremendously. It should not, however, be classed as a signal unless the process is trying to listen for underwater background sound.

As another example, suppose that an attempt is being made to discover the existence of dolphins, and that the boat the hydrophone is mounted on uses electrical power — there will be a high energy band at the 60 Hz line (50 Hz for ships from the USA and some other parts of the world) whenever the ship is under power. Again, this should be classed as noise even though it would traditionally be thought of as a signal against the ocean background.

Finally, consider the same ship discovering a herd of seals. The method should potentially be such that the seals can be ignored on this first occasion, since dolphins are sought, although found again on subsequent trips when seals are of primary interest — luckily this problem can be overcome by separating the problem into identification, discussed in Chapter 6, and classification, discussed briefly in Chapter 7.

These cases illustrate that there are times when signals should be treated as noise for the procedure to be of practical use, or alternatively, sounds should only be classed as signals if they are similar to the object that is being sought.

One approach is to examine the sound at different frequencies, hence the great usage made of Fourier and wavelet techniques in sound analysis, as described in Chapter 2. Examining some different signals, as shown in Figure 3.1, shows that each has differing frequency characteristics when viewed either as a wavelet decomposition or as a spectrogram.

Hence, frequency-based analysis is a good approach to pursue if not all loud events are to be classed as signal. With this in mind, the rest of this section is given over to descriptions of the four classes of events that need to be identified — dolphins, shrimps, seals and rain.

3.3.2 Dolphin Signals

Dolphin “speech” has been studied for several years now (Auty (1986) contains a summary of work by Evans and Prescott in 1962) and three broad classes of sound have been defined based on duration and frequency bands. These are:

(i) Click Duration 0.001 – 0.01 s

Frequency 0.1 – 35 kHz

(ii) Whistle Duration 0.1 – 3.6 s

Frequency 4 – 20 kHz

(iii) Bark  Duration 0.1 s

Frequency 0.2 – 16 kHz

Obviously, several of these sounds can be produced at any time point, and in any combination. Fortunately, the recordings used for this work have several stretches of time containing only clicks or whistles.

Dolphin clicks and whistles are shown in Figure 3.2 as wavelet decompositions and spectrograms. As can be seen, there are visual differences between these sounds when viewed in the frequency domain.

For comparison, Figure 3.3 shows the raw data. The spectrograms are based on 128 frequency bins (see Chapter 2), or 256 time points whilst the wavelet decompositions are based on Daubechies’ order 10 wavelets.

(a) Dolphin whistle  (b) Electric current

(c) Dolphin whistle  (d) Electric current

(e) Dolphin whistle  (f) Electric current

Figure 3.1: Raw signals (top), wavelet decompositions (middle) and spectrograms (bottom) of a dolphin whistle and a 60 Hz electric current.

(a) Dolphin click  (b) Dolphin click

(c) Dolphin whistle  (d) Dolphin whistle

Figure 3.2: Clicks (top) and whistles (bottom) shown as wavelet decompositions (left) and spectrograms (right).

(a) Dolphin click (b) Dolphin whistle

Figure 3.3: Raw plots of dolphin sounds.

One feature of note in the dolphin recordings is the narrow band, high energy signal that is especially apparent on the spectrograms at around the 7 kHz level. This band, and some fainter harmonics in higher frequencies, is due to resonance caused by the recording process. The hydrophone was suspended in the water by a metal rod which started resonating, this being captured by the hydrophone.

As a final word on the dolphin sounds it should be noted that not all species of dolphins make exactly these types of sounds but most fit the broad profile offered, and the Atlantic Bottlenose dolphins presented here are a very good match.

3.3.3 Shrimp Signals

Shrimps generally have less obvious “speech” categories than dolphins, indeed the only audible sound they have been heard to exhibit is a snapping sound that bears many resemblances to the dolphin click, as it is also broad band and of very short duration. Indeed, it takes a highly practised listener to differentiate the two sounds audibly.

Figure 3.4 shows a sequence of shrimp snaps represented in the time domain and two frequency domain representations, the wavelet decomposition and the spectrogram.

The spectrogram is based on 256 data points to give 128 frequency bins whilst the wavelet used is a Daubechies’ order 10. It should also be noted that the 7 kHz resonance line exhibited by the dolphin recording is absent for the shrimp clicks. In this case the hydrophone was suspended by a wire which removes the potential for resonance.

It is actually extremely difficult to find many differences between shrimp snaps and dolphin clicks. As well as being audibly similar, they also exhibit similar properties in the frequency domain (compare Figures 3.2 and 3.4). On exceptionally clear signals however, certain differences can be seen. Figure 3.5 shows a dolphin click and a shrimp snap as raw signals and as wavelet decompositions.

As can be seen, there are some small differences, in particular the shrimp snap behaves far more like a decaying sine function than does the dolphin click, as noted by Smith (1993). It should be remembered however that these examples are atypical in their clarity.

3.3.4 Seal Signals

Seals are in some respects similar to dolphins in that there main characteristic noise, a bark, is similar to the dolphin signal of the same name. As this is the main communication method between seals it does tend to be of a greater amplitude than the dolphin bark, being more akin to the dolphin whistle in this respect.

Figure 3.6 shows a seal bark represented in the time domain and two frequency domain representations, the wavelet decomposition and the spectrogram.

The spectrogram is based on 256 data points to give 128 frequency bins whilst the wavelet used is a Daubechies’ order 10. It should also be noted that the 7 kHz resonance line exhibited by the dolphin recording is also absent for the seal signal In this case the hydrophone was suspended by a wire which removes the potential for resonance.

3.3.5 Falling Rain

The recording gives the sound of rain falling on the surface several metres above the hydrophone in a continuous and essentially random pattern. There are no signals of interest on the recording which makes it ideal as a test data set — any identification method should classify the entire recording as noise.

Figure 3.7 shows a section of rain as raw signal, a wavelet decomposition (based on a Daubechies’ order 10 wavelet) and a spectrogram (based on 128 frequency bins).

Again the recording hydrophone was suspended by wire so there is no resonance band in the 7 kHz range.

(a) The raw signal (b) The wavelet decomposition

(c) The spectrogram

Figure 3.4: Various representations of shrimp snaps.

(a) Shrimp snap  (b) Dolphin click

(c) Shrimp snap  (d) Dolphin click

Figure 3.5: Raw signals (top) and wavelet decompositions (bottom) of a shrimp snap (left) and a dolphin click (right).

(a) The raw signal (b) The wavelet decomposition

(c) The spectrogram

Figure 3.6: Various representations of a seal bark.

(a) The raw signal (b) The wavelet decomposition

(c) The spectrogram

Figure 3.7: Various representations of rain falling.

3.4 The Features of the Artificial Sound

In order to measure how well the identifier performs, and to enable comparisons to be made with other methods, a data set with known signal locations is required.

The recordings displayed above cannot be used for this task as the actual positions of any signals are not known and any method used to identify positions will be subjective, such as audibly identifying signal positions. For this reason, it was decided to construct a data set with similar properties to underwater sound but with known signals at specified positions.

This also allows a range of noise levels (NLs) to be used so that the test can be made gradually more difficult. As will be recalled from Chapter 2, the NL was shown to be given by

NL = 10 log_{10}(\frac{1}{σ^{2}} )

where σ^{2} is the variance of the power of the noise. It was decided to use a variety of NLs from 5 dB down to -10 dB. The noise power was taken to be negative exponential which, from the theoretical results above, should provide a close simulation of underwater sound — this is known as pink noise to differentiate it from white noise in which the power follows a normal distribution.

A signal was then added to this noise. In order to make a fair test the signal needs to be noiseless so that it doesn’t distort the background. It was decided, as all of the underwater recordings were noisy, to use a recording of somebody enunciating the word “operator” under laboratory conditions. Human speech such as this tends to be most similar to dolphin whistles, in that it is narrow band but of long duration.

As human speech is at a much lower frequency band than underwater sound, the recording rate was only 8 kHz, leading to a maximum usable frequency of 4 kHz. This is ideal for human speech although this difference should be noted as it means that the underwater sounds and the operator data are not on comparable scales.

The noise lasts for 4.096 second, i.e. 32768 data points, with the word “operator” being added at the end of this period. This is shown in the following section.

3.4.1 Descriptive Analysis

Figure 3.8 shows the entire length of the operator data at 5 dB as a raw signal, a wavelet decomposition and a spectrogram. In all the subsequent figures, the wavelet decomposition is based around a Daubechies’ Least Asymmetric wavelet of order 8, whilst the spectrogram is based on 128 frequency bins. This choice of wavelet was made after the experiments performed in chapter 4 found it to be sensible choice for many categories of sound.

This is a very clear signal, as it should be for a positive signal to noise ratio, but as the NL is decreased the signal becomes less easy to discern. Indeed, by the time it reaches -8 dB, the signal cannot be audibly distinguished against the background noise. Figures 3.9 and 3.10 show the result of gradually decreasing the noise level.

It should also be noted that the same data set exists with noise levels of 1 dB and −1 dB, although these are not shown here, being visually similar to the 0 dB data set shown.

Obviously, as the NL drops, the actual word “operator” is much harder to visually distinguish. However, the position is still known. Figure 3.11 shows the extreme signal to noise ratios of 5 dB and −10 dB with the actual position of the sound.

3.4.2 Usage of the Operator Data

As the actual position of the signal is known, this data set is ideal for calibrating and testing signal identification methods and for comparing new methods against standard techniques. Indeed, these two tasks are the sole reasons for using this data set.

The ideal position would obviously be if real underwater sound data sets could be used for this task, but this would require knowledge of the exact positions of all the signals in the recording, even those that are visually and

(a) The raw signal (b) The wavelet decomposition

(c) The spectrogram

Figure 3.8: Various representations of the operator data at 5 dB NL.

(a) Operator data with 0 dB NL (b) Operator data with −5 dB NL

(c) Operator data with 0 dB NL (d) Operator data with −5 dB NL

(e) Operator data with 0 dB NL (f) Operator data with −5 dB NL

Figure 3.9: Operator data at 0 and −5 dB NL shown as raw signals (top), wavelet decompositions (middle) and spectrograms (bottom).

(a) Operator data with −8 dB SNR (b) Operator data with −10 dB SNR

(c) Operator data with −8 dB SNR (d) Operator data with −10 dB SNR

(e) Operator data with −8 dB SNR (f) Operator data with −10 dB SNR

Figure 3.10: Operator data at −8 and −10 dB SNR shown as raw signals (top), wavelet decompositions (middle) and spectrograms (bottom).

(a) Operator data at 5 dB SNR

(b) Operator data at −10 dB SNR

(c) Positions of the word operator in the data

Figure 3.11: The position of the word operator in the signal. audibly indistinguishable from the background. Such “marking up” was not possible with the real underwater sound data, hence the necessity of using an artificial data set.

Attempts have been made to make the artificial set as similar as possible to the real situation, hence the use of pink noise for the background and a noiseless signal. However, it should be noted that it is not possible to duplicate the real case exactly and hence there is the potential for differences to exist between results obtained from the artificial and the real situations. These differences should be minimal however.

Now that the data has been described, subsequent chapters turn their attention to techniques for data reduction and signal identification in respect of such data.

Data Compression And Noise Reduction

4.1 Introduction

Data reduction is one of the major aims of signal processing, and is also highly important for the present task. Many signals such as images and the underwater data used in this work, contain large amounts of data which need to be reduced for storage reasons or to make further processing easier. There are many ways in which such data compression can be achieved, although there are two main classes of methods — reduction can be either lossless or lossey. As an example, for compression of image data there are two main algorithms, the first is the JPEG algorithm, mentioned in Chapter 2, that achieves the best compression but is lossey, i.e. it is impossible to recreate the original data exactly from the compressed data as some information is lost, whilst the PNG (Portable Network Graphic) standard achieves less compression but is lossless, i.e. the original data can be obtained from the compressed version.

In the context discussed here, lossey techniques are of primary interest because they remove only unimportant data, retaining significant data which corresponds to the major objective. Wavelet decomposition lends itself to use in lossey compression and will form the focus of the discussion here.

Lossless compression techniques also contain many patent problems whilst lossey techniques tend not to, readers interested in the best lossless techniques should pursue Boutell (1996) which details the PNG compression method. This is considered to be the best available (patent free) method and possibly the best of any method.

The Fourier based JPEG format is the only widely used lossey compression method for image data, although the FBI fingerprint database (see Chapter 2) uses a wavelet based method.

For audio data most lossey compression is based on the Fourier transform, mainly because it is easy to implement and is fast. Less work has been performed on using the wavelet transform as a basis for a compression method although some techniques have been described in the last few years for reducing the amount of non-significant information in such a decomposition.

In this chapter, a wavelet based method has been chosen as the data reduction method for three reasons. Firstly because it will also subsequently be used to identify signals (see Chapter 6) and there is an obvious time saving due to only having to decompose the sound into the wavelet representation rather than into the Fourier domain for the compression and the wavelet domain for signal detection. Second, because the existing wavelet compression methods do not perform well with underwater sounds and a new technique will be presented here that overcomes these. Finally on the basis of the reasons presented in Chapter 2 for using wavelet rather than Fourier techniques in general.

So, the remainder of this chapter will consist of a description of the existing wavelet based methods, a description of the new method, and a comparison between these techniques to show their relative performance.

4.2 Data reduction using wavelets

Wavelet based compression techniques are, by definition, lossey. The wavelet decomposition does not, of itself, perform any data reduction at all. There are exactly as many wavelet coefficients as there were original data points and so data compression can only be achieved by discarding coefficients. How one chooses coefficients to discard is therefore the basis of wavelet compression methods

In any underwater sound, some of the coefficients will correspond wholly or partially to noise rather than signal. As it is necessary to discard some coefficients, it is obviously better to discard those coefficients that contain predominately noise and leave the ones containing signal. Such a technique is known as thresholding and has the advantage that reconstructed sound can not only be compressed but is also less noisy than the original. This is the reason for the primary interest in using a wavelet based compression technique here — subsequent analysis will obviously be easier if there is less noise in the sound.

The next section will detail the general theory behind thresholding of wavelet decompositions and will then provide descriptions of three of the most prominent thresholding methods in use, VisuShrink, SureShrink and Cross Validation as well as briefly mentioning other possible methods before finally presenting RunsThresh, a new addition to these standard thresholding techniques.

4.2.1 Threshold Methods

Donoho and Johnstone (1994) and Donoho et al. (1995) discuss the issue of appropriate thresholds to use for recovering a function of unknown smoothness from noisy, sampled data. If data y_{1}, y_{2}, …, y_{n}(n = 2^{J+1}) can be modelled by

y = f(x_{i}) + ϵ_{i}, i = 1, 2, …, n

where the {ϵ_{i}}_{i=1}^{n} are independent and identically distributed N(0, σ^{2}) and f(·) is the function we would like to estimate, and we use the notation f = {f(x_{i})}_{i=1}^{n} and f̂ = {f̂(x_{i})}_{i=1}^{n} for the vectors of the true and estimated values, then we wish to minimise:

R(f̂, f) = \frac{E∥f̂ − f∥^{2}}{n}

Now, if we let W denote the orthogonal DWT then the white noise above is transformed into white noise (Donoho and Johnstone 1995). Hence the above model is equivalent to

y_{j,k} = w_{j,k} + ϵ_{j,k}, j = 0, 1, …, J; k = 0, 1, …, 2^{j} − 1,

where {y_{j,k}} are the wavelet coefficients of {y_{i}}, {w_{j,k}} are the wavelet coefficients of {f(x_{i})} and {ϵ_{j,k}} is a noise sequence of independent, identically distributed N(0, σ^{2}) variables.

Thus, if ŵ_{j,k} are estimates of the coefficients, then there exists an estimate, f̂ of f which is obtained by f̂ = W^{−1}ŵ, where W^{−1} is the inverse discrete wavelet transform. Therefore, by thresholding the wavelet coefficients w_{j,k} to produce w_{j,k}^{∗}, we can reconstruct an estimate of f by f̂^{∗} = W^{−1}ŵ^{∗}.

Degree of Thresholding

Wavelet based thresholding is based upon the hypothesis that a good wavelet expansion of the data will contain many small coefficients that describe the noise and a few large coefficients that describe the desired signal. From this follows the idea that if we can discover which coefficients correspond to signal, we can set the rest to zero and achieve data reduction and noise reduction by the same process.

Donoho and Johnstone (1994) proposed two degrees of thresholding (i.e. some raw control over the amount of thresholding that should be performed). The most extreme is hard thresholding, whilst the less severe is soft thresholding. These thresholding rules are defined by:

w_{j,k}^{hard} = w_{j,k}I(|w_{j,k}| > λ)
w_{j,k}^{soft} = sign(w_{j,k}) max(0, |w_{j,k}| − λ)

where I is the usual indicator function and λ ≥ 0 is the threshold parameter.

It can be seen that hard thresholding is a straight ‘keep’ or ‘kill’ rule, whilst soft thresholding is a ‘shrink’ or ‘kill’ rule. Hard thresholding is generally better at retaining peak heights and discontinuities whilst soft thresholding gives a smoother result.

From the above expressions, it is readily seen that the value of λ is highly influential. A value that is too extreme will either overfit or underfit the data, producing either too little reduction or a non-recognisable reproduction of the sound. Various methods to estimate λ have been proposed, and several of the more important methods are examined below. First however, a note on the estimation of the noise level, σ, in the data.

The choice of threshold value will depend on this parameter and, as it is highly unlikely to be known, σ must be estimated from the data. The standard deviation of the w_{j,k} provides such an estimate but as Donoho and Johnstone (1994) and Donoho et al. (1995) point out, it is important that the estimate be robust. They propose the use of the median absolute deviation of the wavelet coefficients at the level of the decomposition corresponding to the maximum dilation, divided by 0.6745 (found through repeated testing on various types of sound), as the coefficients at this level are, with few exceptions, pure noise. Hence, we can use the estimate of σ given by

σ̂ = \frac{MAD(w_{J−1,k})}{0.6745}

where J − 1 is the finest level of the decomposition.

The various thresholds will now be examined. Some comparative performance results will be given later but for a simple view of the effects of various threshold methods, the following dolphin whistle (Figure 4.1) will be used.

Figure 4.2 shows a wavelet decomposition of the whistle using a Daubechies’ least asymmetric wavelet of order 8. This is one of the two wavelets used for the majority of the comparative thresholding tests undertaken in this chapter. The other is a Daubechies’ extremal phase wavelet with of order 2.

The descriptions of the various threshold methods will be accompanied, where possible, by comparative diagrams showing the effect of applying soft and hard threshold to the wavelet decomposition of the whistle. More detailed analysis and reconstructions are given later in the comparisons section.

Figure 4.1: A dolphin whistle of length 0.8s.

VisuShrink

Donoho and Johnstone (1994) proposed the universal threshold, which they named VisuShrink. Representing the threshold value by λ_{Visu}, then

λ_{Visu} = σ̂\sqrt{2 log n}

where σ̂ is estimated by the method detailed above and n = 2^{J} (with J being the total number of levels in the decomposition).

Figure 4.3 shows the reconstruction after thresholding by VisuShrink of the whistle in Figure 4.1.

Donoho and Johnstone showed that this estimator is spatially adaptive (it adapts depending on the spread of the data) and is also asymptotically near-minimax over the whole range of Besov spaces. Donoho and Johnstone (1997) have also shown that it asymptotically outperforms all linear estimators (such as splines, kernel estimates and truncated Fourier series) over a range of Besov spaces containing spatially inhomogeneous functions. It is therefore a powerful estimator despite its simplicity. This simplicity does

Figure 4.2: Wavelet decomposition of the whistle.

lead to a further, less desirable, feature however: the estimator is highly conservative as the data is only represented by σ, whatever the form of the function. VisuShrink is in fact conservative enough to asymptotically guarantee that the hard thresholded reconstruction will contain no noise — to achieve this there is also the possibility that some of the signal will also have been removed.

In fact, Donoho and Johnstone suggest that the threshold should only be applied at levels j ≥ j_{0}, where j_{0} is 5, say. This allows the lower frequency bands to be retained even if the coefficients are smaller than the threshold value and goes some way to compensating for the conservative nature of the threshold. A further technique is to use a non-integer resolution level of the wavelet transform, as suggested by Hall and Nason (1997).

This threshold is also near optimal only under the assumption that the noise is i.i.d. N(0, σ^{2}) which, for underwater sounds would be unlikely to be the case. This means that, as the comparisons presented later demonstrate, VisuShrink performs relatively poorly with sounds such as those used in this work.

(a) Hard thresholding  (b) Soft thresholding

Figure 4.3: VisuShrink applied to the whistle.

In an attempt to move away from a such a strict noise model and in a move to formulate less conservative thresholds, several new approaches have been developed. These are all, to some extent, data-driven thresholding rules, some of which will now be discussed.

SureShrink

Donoho and Johnstone (1995) proposed an alternative threshold method, SureShrink, which can be thought of as optimally smoothness-adaptive in that, essentially, a function with jumps will have the jumps retained in the reconstruction and one with smooth sections will also have them retained. This suggests that the advantages of this method are particularly marked when the underlying function has jump discontinuities against a smooth background.

Figure 4.4 shows the reconstruction of the whistle after thresholding by SureShrink (which can only use the soft thresholding variety). The actual threshold is described below but it should immediately be noted that the SureShrink threshold is level-dependent, that is, the value of λ is different for each level, j, of the decomposition, λ_{j} say.

Figure 4.4: SureShrink applied to the whistle.

The actual threshold, λ_{Sure_{j}} is given by

λ_{Sure_{j}} = argmin_{0≤λ≤\sqrt{2 log(n_{j})}}SURE(λ_{j}; x)

where n is the number of data points, j is the resolution level and SURE(λ; x) is given in the Donoho and Johnstone (1995) but is essentially Stein’s unbiased risk estimator, after Stein (1981). In the paper, it is shown that this is a better threshold for real data than VisuShrink in many situations. The calculation of this threshold is an O(n log(n)) process which is close to that for VisuShrink.

The actual SureShrink estimator is slightly more complicated in practice as extremely sparse wavelet coefficients lead to a distorted SURE profile. This leads to the use of a hybrid scheme which overcomes this problem. Such a scheme, as shown in Donoho and Johnstone (1995) uses a VisuShrink threshold when the data is very sparse and the SureShrink threshold otherwise.

SureShrink is also shown to be asymptotically near-minimax and a practical demonstration of the advantages of using this threshold approach is given by Johnstone and Silverman (1997).

Cross Validation

Nason (1996) has suggested that a cross-validation approach should be adapted to choose the threshold. This work has been extended by Wang (1996) and Weyrich and Warhola (1994).

Classical cross validation relies on some form of ‘leave one out’ scheme, where one data item is dropped and some statistic calculated over the remaining data items. Obviously, under a wavelet decomposition, such a scheme is impractical as the fast wavelet transform relies on having a power of two data items. The solution is to drop half the data items when cross validating to go from 2^{j} to 2^{j−1} data items. This is called two-fold cross-validation and is described in Nason (1996) and Nason (1997). The latter work also has descriptions of the extensions implemented by Wang (1996) and Weyrich and Warhola (1994).

The aim of this threshold is to minimise the MISE (Mean Integrated Square Error) between the true function and the thresholded wavelet estimator. That is, the threshold t should minimise

M(t) = E ∫ {f̂(x) − f(x)}^{2}dx.

Obviously, the true function is not usually known when thresholding, so an estimate of M must be devised. This estimate is formed from the two sets of data items obtained by splitting the original data set in two.

From this estimate, a cross validated threshold value can be obtained. This process is O(n) and simulation studies have shown it to be a good threshold method (Abramovich et al. 1998).

Figure 4.5 shows the reconstructions of the whistle after thresholding by Cross Validation.

Figure 4.5: Cross Validation applied to the whistle.

Other Methods

The search for a better thresholding method has produced various other developments.

One is the approach taken by Vidakovic (1994), Clyde et al. (1995) and Chipman et al. (1996), and shown in more detail by Abramovich et al. (1998), where a Bayesian framework is considered. The prior is considered as a mixture of normal values with a magnitude of the large wavelet coefficients, and point masses at zero.

A further approach has been to consider thresholding as a multiplehypothesis testing procedure where blocks of coefficients are tested to see if they equal zero simultaneously. This technique has been developed by Abramovich and Benjamini in their 1995 and 1996 papers (where they called it the False Discovery Rate estimator) and by Ogden and Parzen (1996b, 1996a).

Recently, Johnstone and Silverman (1997) have developed a technique for data with correlated noise based on a level-dependent threshold methods. Some of the above approaches can be extended using these techniques to cope with such noise structures.

4.2.2 RunsThresh

One of the first areas considered at the start of this work was that of data reduction (with sound enhancement as a by-product of such techniques). At the time, the best known thresholds were VisuShrink and SureShrink. After investigation of these methods, it was felt that neither were suitable for the particular data used in this problem, as will be illustrated later in this chapter.

To overcome the deficiencies exposed in these standard methods, a modification of VisuShrink was proposed. This modification is based on the standard non-parametric runs-test, described, for example, by Siegel (1956) and hence is known as RunsThresh. It is briefly described in Powell et al. (1995) but will be discussed in more detail here.

If a section of the wavelet decomposition in Figure 4.2 is enlarged, as shown in Figure 4.6, then a sinusoidal pattern in the coefficients is clearly visible.

Figure 4.6: Magnified section of the whistle wavelet decomposition.

This feature is the basis of the reasoning for using the runs test. Most data do not exhibit this pattern in the wavelet coefficients, making it reasonable to assume that they are, in some way, related to the unique features of the underwater signals. This being the case, any threshold should ideally retain these coefficients when they are present. This can be thought of as an empirical way to move away from the conservative nature of the VisuShrink threshold.

The runs test was considered the most appropriate technique to use in discovering if such patterns exist in a level of wavelet coefficients (making the technique, by necessity, level-dependent) and the following scheme was proposed as an adaption of the VisuShrink technique.

For each level of the wavelet decomposition, the coefficients are grouped into a set time period, which is identical for all levels. This means that levels corresponding to higher frequencies have more coefficients in a group. If a level has groups of less than 32 coefficients then the VisuShrink threshold itself is used (this is a similar idea to that of only applying the VisuShrink threshold for lower levels of the decomposition, as described earlier).

Now, for all groups of 32 or more coefficients, the significance of runs within each group is assessed on the basis of the observed number of runs of positive and negative coefficients, by way of the usual approximate test statistic, as in Siegel (1956):

z = \frac{r − (\frac{2n_{1}n_{2}}{n_{1}+n_{2}} + 1)}{\sqrt{\frac{2n_{1}n_{2}(2n_{1}n_{2}−n_{1}−n_{2})}{(n_{1}+n_{2})^{2}(n_{1}+n_{2}−1)}}}

where r is the number of runs, n_{1} is the number of positive values, n_{2} the number of negative values and z has an approximately N(0, 1) distribution if the runs are randomly distributed.

If the two-tailed test is significant at the 5% level, then the threshold used, λ_{Runs} is given by

λ_{Runs} = κλ_{Visu}

where λ_{Visu} is the VisuShrink threshold as defined above and 0 < κ < 1 is a scaling factor. If the test is not significant then κ is set to 1.

Now, as can be seen from the above explanation, several parameters can effectively be used to control the amount of thresholding performed. The size of a coefficient block to runs test, κ and the significance level tested at will all determine the threshold level. Chapter 5 examines these parameters in more depth, as well as providing more evidence for the effectiveness of RunsThresh with underwater sounds.

For now, it will be noted that the values used to produce the results in the remainder of this work are a block of coefficients with length 0.00625 seconds, a κ value of 0.1 and a significance level of 5%.

It is also, of course, possible to apply a similar technique in conjunction with any other threshold method, such as SureShrink by replacing λ_{Visu} with λ_{Sure} in the above equations.

Therefore, if there is a pattern of runs (such as where there is a sinusoidal pattern in the coefficients), the threshold is set to a lower value than is normal and more of the coefficients are retained, giving the desired effect.

The reconstructions of the whistle after thresholding by the VisuShrink based RunsThresh threshold are shown in Figure 4.1 are given in Figure 4.7. As can be seen, the prominent sinusoidal patterns have been retained whilst other features have either been shrunk or killed.

(a) Hard thresholding  (b) Soft thresholding

Figure 4.7: RunsThresh applied to the whistle.

This threshold policy has been implemented as a new option for the threshold.wd() function in the WaveThresh package for S-Plus (Nason 1993), as described in Nason and Silverman (1994). The syntax is fully detailed in appendix A but simply consists of three additional flags.

4.3 Comparisons

The differing abilities of these thresholds will now be considered. In common with other papers on thresholding, we will look at the four functions defined in, for example, Donoho and Johnstone (1995). These are Blocks, Bumps, HeaviSine and Doppler. The various functions are defined in Table 4.1.

Blocks

f(t) = ∑ h_{j}K(t − t_{j}) K(t) = (1 + sign(t))/2
(t_{j}) = (.1, .13, .15, .23, .25, .40, .44, .65, .76, .78, .81)
(h_{j}) = (4, −5, 3, −4, 5, −4.2, 2.1, 4.3, −3.1, 2.1, −4.2)

Bumps

f(t) = ∑ h_{j}K(\frac{t−t_{j}}{w_{j}}) K(t) = (1 + |t|)^{−4}
(t_{j}) = t_{Blocks}
(h_{j}) = (4, 5, 3, 4, 5, 4.2, 2.1, 4.3, 3.1, 5.1, 4.2)
(w_{j}) = (.005, .005, .006, .01, .01, .03, .01, .01, .005, .008, .005)

HeaviSine

f(t) = 4 sin 4πt − sign(t − .3) − sign(.72 − t)

Doppler

f(t) = \sqrt{t(1 − t)} sin(2π \frac{1−ϵ}{t+ϵ}), ϵ = .05

Table 4.1: Formulae for test functions

These functions are shown in Figure 4.8; all have 1024 data points. It should be noted that these figures, while representing some real life phenomena, do not accurately represent underwater signals. Their purpose here is to provide some common reference points for comparison of RunsThresh with other thresholds.

As RunsThresh is designed specifically for underwater sound, a more realistic test of the quality of the thresholding will also be given by using such sounds. For this section, a selection of typical short signals have been extracted. These are shown in Figure 4.9.

In common with Abramovich et al. (1998), 100 replicates of the test functions were generated by adding independent random noise, ϵ_{i} ∼ N(0, σ^{2})

(a) Blocks  (b) Bumps

(c) HeaviSine  (d) Doppler

Figure 4.8: The four test functions defined in Table 4.1

to the 1024 data points on the interval [0, 1]. The values of σ were taken so that the noise levels were SD(f)/σ = 10, 7, 5 and 3. Similarly, the wavelets used were Daubechies’ least asymmetric wavelets of order 8 and her extremal phase of order 2. Hard and soft forms of the thresholds were compared using a goodness-of-fit measured by the average mean square error (AMSE):

AMSE(f̂; f) = \frac{1}{n}E∥f̂− f∥^{2} = \frac{1}{n} ∑_{i=1}^{n}(f̂_{i} − f_{i})^{2}

(a) Dolphin Whistle  (b) Seal Bark

(c) Shrimp Clicks  (d) Rain

Figure 4.9: The four typical sounds used for the comparisons

Figure 4.10 shows the four test functions with added noise of level 3, this is the noisiest of the four possible levels and all future figures will also use noise level 3.

The full results of these analyses are given in Tables 4.2 through 4.5. For comparison with Abramovich et al. (1998), all tables show the AMSEs for the 100 replicates using either hard or soft thresholding and Daubechies’ Least Asymmetric wavelets of order 8 or her Extremal Phase wavelets of order 2.

Two forms of RunsThresh results are presented, the first RunsVisu is based on VisuShrink whilst the second, RunsSure is based on SureShrink.

(c) HeaviSine  (d) Doppler

Figure 4.10: The four test functions with added noise

The full results tend to indicate that soft cross validation is the best method for these test signals, although is should be noted that for many noise/signal combinations, hard RunsVisu thresholding is on a par with it.

It should of course be noted that in general, the runs and non-runs thresholded versions tend to give very similar results. This is not overly surprising as a difference should only be discernible if the signal has significant runs — this is not generally the case in these test signals.

It is also helpful to view the various reconstructions after thresholding.

SNR Blocks  Bumps  HeaviSine Doppler
RunsVisu  10 1.31 (0.113) 1.71 (0.187) 0.11 (0.015) 0.46 (0.054)
7 2.12 (0.160) 3.06 (0.326) 0.17 (0.028) 0.78 (0.086)
5 3.22 (0.207) 4.97 (0.443) 0.23 (0.050) 1.22 (0.148)
3 5.26 (0.312) 10.15 (0.817) 0.43 (0.145) 2.38 (0.233)
VisuShrink  10 1.36 (0.099) 1.86 (0.149) 0.10 (0.009) 0.47 (0.046)
7 2.24 (0.153) 3.28 (0.270) 0.15 (0.016) 0.78 (0.066)
5 3.27 (0.205) 5.23 (0.416) 0.20 (0.018) 1.22 (0.111)
3 5.27 (0.284) 10.40 (0.689) 0.32 (0.043) 2.33 (0.189)
SureShrink  10 0.38 (0.047) 0.43 (0.055) 0.11 (0.009) 0.16 (0.024)
7 0.56 (0.084) 0.67 (0.108) 0.16 (0.015) 0.32 (0.210)
5 2.59 (1.607) 0.85 (0.062) 0.21 (0.017) 1.38 (0.102)
3 5.84 (0.268) 2.66 (0.302) 0.32 (0.043) 2.52 (0.167)
Cross  10 0.23 (0.019) 0.25 (0.020) 0.06 (0.007) 0.11 (0.012)
Validation  7 0.41 (0.032) 0.49 (0.032) 0.10 (0.014) 0.21 (0.023)
5 0.72 (0.057) 0.85 (0.053) 0.23 (0.027) 0.38 (0.044)
3 1.66 (0.129) 2.07 (0.149) 0.31 (0.049) 0.89 (0.096)
RunsSure  10 0.38 (0.050) 0.43 (0.054) 0.12 (0.014) 0.17 (0.027)
7 0.58 (0.085) 0.68 (0.112) 0.17 (0.026) 0.34 (0.206)
5 2.56 (1.557) 0.94 (0.108) 0.23 (0.047) 1.37 (0.146)
3 5.81 (0.300) 2.83 (0.367) 0.43 (0.139) 2.55 (0.217)

Table 4.2: AMSE for the different threshold types using standard test signals, with standard errors in parentheses. Daubechies’ Least Asymmetric wavelets of order 8 and soft thresholding used.

SNR Blocks  Bumps  HeaviSine Doppler
RunsVisu  10 1.00 (0.105) 1.80 (0.180) 0.12 (0.014) 0.85 (0.068)
7 1.74 (0.160) 3.10 (0.322) 0.18 (0.029) 1.31 (0.091)
5 2.78 (0.226) 4.94 (0.447) 0.26 (0.061) 1.96 (0.126)
3 5.07 (0.404) 9.41 (0.805) 0.43 (0.127) 3.45 (0.279)
VisuShrink  10 1.02 (0.100) 1.84 (0.172) 0.11 (0.009) 0.84 (0.067)
7 1.76 (0.155) 3.18 (0.299) 0.16 (0.018) 1.30 (0.084)
5 2.80 (0.217) 5.07 (0.364) 0.23 (0.019) 1.94 (0.124)
3 5.05 (0.389) 9.63 (0.720) 0.35 (0.044) 3.34 (0.243)
SureShrink  10 0.26 (0.034) 0.39 (0.058) 0.12 (0.009) 0.29 (0.033)
7 0.44 (0.067) 0.61 (0.100) 0.17 (0.018) 0.49 (0.280)
5 1.51 (1.349) 0.82 (0.069) 0.23 (0.018) 2.21 (0.116)
3 5.64 (0.323) 2.98 (1.735) 0.35 (0.043) 3.64 (0.204)
Cross  10 0.21 (0.025) 0.37 (0.053) 0.06 (0.008) 0.18 (0.011)
Validation  7 0.34 (0.031) 0.54 (0.052) 0.10 (0.015) 0.33 (0.022)
5 0.57 (0.048) 0.86 (0.079) 0.16 (0.024) 0.58 (0.042)
3 1.42 (0.111) 2.00 (0.176) 0.33 (0.049) 1.26 (0.104)
RunsSure  10 0.27 (0.035) 0.40 (0.059) 0.13 (0.013) 0.30 (0.038)
7 0.45 (0.070) 0.62 (0.101) 0.19 (0.027) 0.52 (0.281)
5 1.52 (1.310) 0.86 (0.096) 0.26 (0.059) 2.23 (0.125)
3 5.65 (0.338) 3.08 (1.740) 0.42 (0.119) 3.74 (0.236)

Table 4.3: AMSE for the different threshold types using standard test signals, with standard errors in parentheses. Daubechies’ Extremal Phase wavelets of order 2 and soft thresholding used.

SNR Blocks  Bumps  HeaviSine Doppler
RunsVisu  10 0.43 (0.047) 0.47 (0.050) 0.10 (0.030) 0.16 (0.033)
7 0.84 (0.124) 0.82 (0.088) 0.14 (0.059) 0.30 (0.068)
5 1.63 (0.223) 1.53 (0.182) 0.27 (0.139) 0.59 (0.122)
3 3.55 (0.438) 4.33 (0.520) 0.64 (0.404) 1.49 (0.393)
VisuShrink  10 0.42 (0.045) 0.47 (0.046) 0.07 (0.007) 0.14 (0.021)
7 0.82 (0.113) 0.81 (0.082) 0.10 (0.018) 0.28 (0.042)
5 1.59 (0.189) 1.52 (0.192) 0.19 (0.030) 0.52 (0.069)
3 3.36 (0.336) 4.19 (0.460) 0.32 (0.043) 1.32 (0.233)
Cross  10 0.23 (0.035) 0.28 (0.089) 0.06 (0.013) 0.10 (0.019)
Validation  7 0.47 (0.112) 0.49 (0.092) 0.11 (0.022) 0.17 (0.027)
5 0.86 (0.197) 0.87 (0.267) 0.16 (0.044) 0.34 (0.066)
3 2.11 (0.351) 2.13 (0.367) 0.36 (0.099) 0.85 (0.167)

Table 4.4: AMSE for the different threshold types using standard test signals, with standard errors in parentheses. Daubechies’ Least Asymmetric wavelets of order 8 and hard thresholding used.

SNR Blocks  Bumps  HeaviSine Doppler
RunsVisu  10 0.22 (0.042) 0.41 (0.041) 0.09 (0.032) 0.36 (0.047)
7 0.56 (0.100) 0.84 (0.102) 0.15 (0.065) 0.59 (0.079)
5 1.20 (0.193) 1.77 (0.214) 0.27 (0.112) 0.94 (0.159)
3 2.72 (0.407) 4.08 (0.395) 0.57 (0.327) 2.30 (0.452)
VisuShrink  10 0.21 (0.034) 0.41 (0.041) 0.08 (0.011) 0.34 (0.035)
7 0.53 (0.093) 0.84 (0.095) 0.11 (0.016) 0.56 (0.062)
5 1.15 (0.155) 1.74 (0.207) 0.19 (0.036) 0.86 (0.080)
3 2.58 (0.321) 3.92 (0.320) 0.35 (0.042) 2.10 (0.260)
Cross  10 0.40 (0.355) 0.66 (0.499) 0.06 (0.010) 0.19 (0.022)
Validation  7 0.50 (0.324) 0.70 (0.446) 0.10 (0.017) 0.35 (0.042)
5 0.63 (0.221) 1.04 (0.391) 0.17 (0.035) 0.61 (0.079)
3 1.48 (0.191) 2.22 (0.403) 0.38 (0.091) 1.38 (0.275)

Table 4.5: AMSE for the different threshold types using standard test signals, with standard errors in parentheses. Daubechies’ Extremal Phase wavelets of order 2 and hard thresholding used.

As the results show, there is a difference between hard and soft thresholding. However, the figures shown here concentrate on soft thresholded reconstructions — the differences with hard thresholding are not important for these visual comparisons. The figures are based around decompositions using the Daubechies Least Asymmetric wavelet of order 8.

Figure 4.11 shows the results for RunsThresh, specifically the VisuShrink based variant. Figure 4.12 shows the results for VisuShrink, Figure 4.13 gives the SureShrink results and finally, Figure 4.14 the Cross-validation results. All of these reconstructions are based on soft thresholding.

(a) Blocks  (b) Bumps

(c) HeaviSine  (d) Doppler

Figure 4.11: Test signals after RunsThresh has been applied

(c) HeaviSine  (d) Doppler

Figure 4.12: Test signals after VisuShrink has been applied

(c) HeaviSine  (d) Doppler

Figure 4.13: Test signals after SureShrink has been applied

(c) HeaviSine  (d) Doppler

Figure 4.14: Test signals after Cross-validation has been applied

Subject 1 Subject 2 Subject 3 Subject 4 Average
RunsVisu  2  1  1  2  1
RunsSure  1  2  3  1  2
Cross-validation  3  3  2  4  3
SureShrink  5  4  4  3  4
VisuShrink  4  5  5  5  5

Table 4.6: Rankings of audible quality after reconstruction of various underwater sounds. Daubechies’ Least Asymmetric wavelets of order 8 and soft thresholding used

Subject 1 Subject 2 Subject 3 Subject 4 Average
RunsVisu  1  2  1  1  1
Cross-validation  2  1  2  2  2
VisuShrink  3  3  3  3  3

Table 4.7: Rankings of audible quality after reconstruction of various underwater sounds. Daubechies’ Least Asymmetric wavelets of order 8 and hard thresholding used

As was noted above, these functions provide the standard test for new threshold methods and so it is no surprise that RunsThresh does not exceed the standards of its analogous standard thresholds. RunsThresh was developed solely for the purpose of thresholding underwater sounds so a more interesting analysis is what happens when all the methods are applied to the test sounds shown earlier.

The problem that occurs when trying to analyse such results however, is that noise already exists in the decomposed section of data, so using differing noise levels is highly problematic. A further, but related, difficulty, is that the true signal is not known so it is impossible to measure how close any reconstruction is to the true (“clean”) signal.

The method proposed by Powell et al. (1995) was to ask a panel of musicians to rank the reconstructions. The original was played before each reconstruction, with each member of the panel being played the reconstructions in a different order to remove any inherent bias. Any of the reconstructions could be replayed on demand and the final rankings are given in Tables 4.6 through 4.9.

Subject 1 Subject 2 Subject 3 Subject 4 Average
RunsVisu  1  2  1  1  1
RunsSure  3  1  3  3  3
Cross-validation  2  3  2  2  2
SureShrink  4  4  5  4  4
VisuShrink  5  5  4  5  5

Table 4.8: Rankings of audible quality after reconstruction of various underwater sounds. Daubechies’ Extremal Phase wavelets of order 2 and soft thresholding used

Subject 1 Subject 2 Subject 3 Subject 4 Average
RunsVisu  1  1  1  1  1
Cross-validation  2  2  3  2  2
VisuShrink  3  3  2  3  3

Table 4.9: Rankings of audible quality after reconstruction of various underwater sounds. Daubechies’ Extremal Phase wavelets of order 2 and hard thresholding used

As these results show, RunsVisu performs considerably better than other methods for underwater sounds and hence in future we will use RunsThresh to refer to a this variety of threshold. These results are confirmed by looking at the reconstructions. Figure 4.15 shows the RunsThresh results, whilst Figures 4.16, 4.17 and 4.18 are for the VisuShrink, SureShrink and Crossvalidation thresholds respectively. Again, all reconstructions are based on soft thresholding.

Given that RunsThresh has been shown to be the best performer for underwater sounds, it is now important to discover if any particular wavelet is best as a basis for the decomposition. A similar test was used as for the choice of threshold but with different wavelets being used for the decompositions and RunsThresh being applied in all cases. Only the least asymmetric and extremal phase wavelets were considered, but orders were varied and both dolphins and shrimp signals were used. The results are presented in full in Table 4.10.

These results would indicate that the choice of wavelet did not make a large impact on the quality of reconstruction, the absolute difference in scores

(c) Shrimp Click  (d) Rain

Figure 4.15: Underwater signals after RunsThresh has been applied

Wavelet Subject 1 Subject 2 Subject 3 Subject 4 Average
Whistle LA 8  1  2  5  6  2
Whistle EP 2  3  1  7  5  3
Whistle EP 6  7  3  6  1  4
Whistle LA 4  2  6  1  4  1
Shrimps LA 8  5  4  2  7  5
Shrimps EP 2  8  8  3  2  6
Shrimps EP 6  4  7  4  8  8
Shrimps LA 4  6  5  8  3  7

Table 4.10: Rankings of audible quality after reconstructing of various underwater sounds using various wavelets and hard thresholding with RunsThresh

(c) Shrimp Click  (d) Rain

Figure 4.16: Underwater signals after VisuShrink has been applied

(c) Shrimp Click  (d) Rain

Figure 4.17: Underwater signals after SureShrink has been applied

(c) Shrimp Click  (d) Rain

Figure 4.18: Underwater signals after Cross-validation has been applied

Threshold  Type Compression Ratio
VisuShrink Hard 74.0%
Soft  74.0%
SureShrink Soft 9.3%
Cross Validation Hard  20.1%
Soft  70.1%
RunsThresh Hard 39.3%
Soft  39.3%

Table 4.11: Compression ratios different thresholds methods

being minimal. So the choice of wavelet is not particularly significant for underwater sounds, or indeed for the test signals shown earlier. Given this observation, the rest of this work will adopt the Daubechies’ Least Asymmetric of order 8 for consistency with other texts on this subject.

Finally, given that one of the initial reasons for thresholding is data reduction, some results for the amount of compression achievable are given. These results were calculated using the percentage of zeros in the wavelet decomposition objects after thresholding. By comparison, the wavelet decomposition of the original whistle had no zeros. Under this system, large compression ratios are best.

Table 4.11 gives compression ratios for various threshold methods using a Daubechies’ Least Asymmetric wavelet of order 8. As can be seen, RunsThresh gives reasonable compression ratios, a product of its noise reduction abilities. Although cross validation and VisuShrink give more impressive compression rates, they do so at the expense of reconstruction quality. As the main aim for this part of the work is noise reduction, with compression being a secondary aim, it was felt to be more useful to achieve reasonable compression rates (such as the 40% achieved by RunsThresh) with high quality noise reduction.

Table 4.12 shows the compression rates achieved using Daubechies’ Least Asymmetric wavelet of order 8 with both hard and soft RunsThresh applied to the standard underwater sounds used throughout this work. As can be seen, the different categories of underwater sound are all compressed to about the same degree.

Signal Type  Hard Thresholding Soft Thresholding
Dolphin Whistle  39.9%  39.3%
Dolphin Clicks  34.0%  34.0%
Shrimp Clicks  65.2%  65.2%
Seal Bark  41.3%  41.3%
Rain  38.3%  36.3%

Table 4.12: Compression ratios for different underwater sounds

RunsThresh - Parameters and Effects

5.1 Introduction

In Chapter 4, a modified thresholding method called RunsThresh was introduced. The results in that chapter demonstrated that RunsThresh has potential in respect of the analysis of underwater sounds. It was also noted that the procedure involves several parameters, the significance level, the block size and κ (the threshold effect) although these were not discussed in any detail.

The purpose of this chapter is to demonstrate that this method does produce significantly better results than less complex adaptations of VisuShrink and to investigate the differences that varying these parameters produces.

Given these aims, we will first investigate the improvements RunsThresh offers and secondly analyse the parameters to discover the best values for our particular underwater sounds.

5.2 The RunsThresh Effect

A possible criticism that could be levelled at the RunsThresh method described in Chapter 4 is that it is simply VisuShrink but with a new, lower, threshold chosen universally. To counter this argument we must show that the lower threshold level is applied selectively, that is, not all blocks of coefficients have significant sets of runs in them.

Consider the whistle used in Chapter 4 to demonstrate various threshold techniques (Figure 4.1). Figure 5.1 shows the blocks of coefficients found to be significant when this is subjected to the hard and soft varieties of RunsThresh and also confirms that there is very little difference between hard and soft thresholding with RunsThresh. If RunsThresh were simply VisuShrink with a lower threshold we would expect all coefficient blocks to be significant.

As can be seen, only some of the blocks are significant and would therefore benefit from the reduced threshold level. In particular, blocks in all levels apart from 12, 13 and 14 are not significant whereas most blocks in these rows are significant. This corresponds well with the wavelet decomposition of the whistle which shows the characteristic features to be concentrated in these three levels.

Figure 5.2 shows the results of thresholding with RunsThresh and with VisuShrink where the VisuShrink threshold has been set to 10% of the normal value and RunsThresh has a κ value of 0.1. The depiction of the wavelet decompositions here differ from those used in the rest of this thesis, in that each scale level has been scaled individually so that the differences between VisuShrink and RunsThresh can be readily discerned.

It should be noted that the figures produced after thresholding show the reduced VisuShrink decomposition contains 33% more non-zero coefficients than the RunsThresh decomposition (30% compression compared with 39%).

RunsThresh does therefore behave in a different manner to VisuShrink. However, the properties of RunsThresh must also be investigated in more detail, in particular to decide on appropriate levels for the different parameters available in our particular application.

The next section examines the effect of these various parameters on the thresholding results.

(a) Hard thresholding (b) Soft thresholding

Figure 5.1: Significant blocks of coefficients (shown in black) during RunsThresh.

(a) RunsThresh (b) VisuShrink

Figure 5.2: The wavelet decomposition of the whistle after thresholds have been applied.

5.3 The RunsThresh Parameters

As noted above, the RunsThresh procedure is controlled by three parameters, the significance level used for the runs-test, the block size and the parameter κ which controls the proportion of the VisuShrink threshold used for a significant block of coefficients.

To analyse the effects of changing these parameters we rely on two techniques. The initial method examines how a set of summary measures alters with the changing parameter value. These summary measures seek to represent some of the features that one would hope to extract from a good threshold. The second method is that used in Chapter 4, namely audio analysis of reconstructions by a panel of musicians. As explained, this method is the best that we could devise to capture the full spectrum of threshold effects.

These techniques were applied to two of the sounds shown in earlier chapters, the dolphin whistle and the shrimp clicks (as these represent the two main classes of signal). The results for these sounds will then be used to choose values for the analysis performed elsewhere in this work and explain why certain parameter values were chosen when a range of values could have been used.

Before the results for the different parameters are presented we will first explain the summary measures used in the comparisons. One of these is also used in Chapter 4 — the compression ratio. This is the amount of the sound that has been set to zero as a percentage of the total sound. This obviously means that high values can be considered good.

The remaining measures are not used elsewhere in this work. They summarise differences between the original sound and the thresholded sound. As explained in Chapter 4, since the values for the original sound with no noise are not known, these measures do not necessarily show if a thresholded version of a sound is good or bad. They can be useful however in showing the effect of changing parameters when used in conjunction with the panel test, where subjective judgements are made about the quality of sound reproduction.

The first measure is the sum of squares of coefficient differences between the thresholded sound and the original, non-thresholded sound. This is given by:

SS_{tot} = ∑_{j} ∑_{k}(d_{j,k}^{∗} − d_{j,k})^{2}

where d_{j,k}^{∗} are the coefficients of the thresholded wavelet decomposition and d_{j,k} are the coefficients of the non-thresholded wavelet decomposition as given in Chapter 2. This measure increases as the number of non-zero coefficients is reduced but, as a raw number, it is difficult to determine whether a particular value is good or not. To help in this task, the second measure is this value expressed as a percentage of the value produced when VisuShrink is used:

%SS_{tot} \frac{∑_{j} ∑_{k}(d_{j,k}^{∗} − d_{j,k})^{2}}{∑_{j} ∑_{k}(d_{j,k}^{v} − d_{j,k})^{2}} ∗ 100

where d_{j,k}^{∗} and d_{j,k} are as given above and d_{j,k}^{v} are the coefficient values for a VisuShrink wavelet decomposition.

The next set of values are the sum of squares of coefficient differences within each scale level. These measures allow the effect of the threshold at certain frequencies to be seen which gives a more balanced picture of the behaviour of the threshold. For scale j, the measure is given by:

SS_{j} = ∑_{k}(d_{j,k}^{∗} − d_{j,k})^{2}

These measures will also increase as more coefficients are reduced, but they will more easily show which parts of the sound are being thresholded. Finally, these measures are also given as a percentage of the VisuShrink values to allow comparisons to be performed more easily. At scale level j, such a parameter would be given by:

%SS_{j} = \frac{∑_{k}(d_{j,k}^{∗} − d_{j,k})^{2}}{∑_{k}(d_{j,k}^{v} − d_{j,k})^{2}} ∗ 100

It should be noted that the values of SS_{j} and %SS_{j} are only given for

κ  0.1  0.3  0.5  0.8
Compression ratio 39.32% 53.64% 62.94% 71.09%
SS_{tot}(×10^{8})  22.97 29.66 46.95 87.23
%SS_{tot}  19.97% 25.78% 40.81% 75.82%
SS_{10}(×10^{8})  2.05  2.05  2.05  2.05
SS_{11}(×10^{8})  9.48  9.48  9.48  9.48
SS_{12}(×10^{8})  4.44  4.68  5.38  9.09
SS_{13}(×10^{8})  0.05  1.18  4.85  15.93
SS_{14}(×10^{8})  0.56  5.89  18.81 44.30
%SS_{10}  100% 100% 100% 100%
%SS_{11}  100% 100% 100% 100%
%SS_{12}  35.69% 37.57% 43.21% 72.97%
%SS_{13}  0.18% 4.11% 16.87% 55.42%
%SS_{14}  1.00% 10.53% 33.63% 79.19%

Table 5.1: Values of summary measures after reconstruction of whistles barks. Daubechies’ Least Asymmetric wavelets of order 8 and soft thresholding used

j = 10 … 14 as in the remaining levels there are too few coefficients for the runs-test to be applied. This means that these levels are simply thresholded with VisuShrink.

Given these summary measures, we will now examine each of the parameters for RunsThresh and determine suitable values to use in the rest of this thesis.

5.3.1 κ, the threshold effect

The threshold effect, κ, runs over the range 0 to 1. A value of 1 produces the universal threshold for all coefficients whilst a value of 0 leaves all significant blocks of coefficients unthresholded.

A suitable value of κ needs to be found such that the amount of thresholding for coefficients in a significant group of runs is small enough that large coefficients from the group will not be removed.

The summary measures described earlier for different values of κ are shown in Tables 5.1 and 5.2.

These results were produced with a block size of 0.00625 seconds and a

κ  0.1  0.3  0.5  0.8
Compression ratio 65.18% 75.8% 80.94% 84.79%
SS_{tot}(×10^{8})  14.76 16.14 18.80 24.03
%SS_{tot}  53.59% 58.6% 68.26% 87.23%
SS_{10}(×10^{8})  1.61  1.61  1.61  1.61
SS_{11}(×10^{8})  3.29  3.29  3.29  3.29
SS_{12}(×10^{8})  5.22  5.22  5.22  5.22
SS_{13}(×10^{8})  0.05  1.04  3.49  8.42
SS_{14}(×10^{8})  2.81  3.20  3.41  3.69
%SS_{10}  100% 100% 100% 100%
%SS_{11}  100% 100% 100% 100%
%SS_{12}  99.4% 99.4% 99.4% 99.55%
%SS_{13}  0.41% 8.76% 29.36% 70.98%
%SS_{14}  75.08% 85.5% 91.26% 98.68%

Table 5.2: Values of summary measures after reconstruction of shrimp clicks. Daubechies’ Least Asymmetric wavelets of order 8 and soft thresholding used

significance level of 5%. This block size limits RunsThresh to levels above 11. Equally, as explained in Chapter 4, levels below 5 are not thresholded at all.

In general, the shrimp click data is much closer to the VisuShrink results than the dolphin whistle data. This is not surprising as the shrimp clicks have far fewer runs of data than the dolphin clicks, indeed, as explained in Chapter 4, RunsThresh was developed from features initially observed in the dolphin whistle wavelet decompositions.

As can be seen, varying κ makes very little difference to level 12 — with the shrimp clicks all the values are close to the VisuShrink level whilst for the dolphin data, only a κ value of 0.8 is significantly different to other κ values. However, for levels 13 and 14, changing the κ value makes a major difference to the results. For the dolphin data, κ values up to 0.5 provide less than 35% of the sum of squared differences that VisuShrink does. For the shrimp clicks, only level 13 provides this level of difference, level 14 results are far closer to the VisuShrink results.

In general, we can see that as κ is increased, the results converge to the VisuShrink results as expected. Indeed, when κ = 0.8, the compression ratio

κ Subject 1 Subject 2 Subject 3 Subject 4 Average
0.8  4  4  4  3  4
0.5  3  3  3  4  3
0.3  2  2  1  2  2
0.1  1  1  2  1  1

Table 5.3: Rankings of audible quality after reconstruction of dolphins whistles. Daubechies’ Least Asymmetric wavelets of order 8 and soft thresholding used

κ Subject 1 Subject 2 Subject 3 Subject 4 Average
0.8  3  4  3  4  4
0.5  4  3  2  3  3
0.3  2  2  4  2  2
0.1  1  1  1  1  1

Table 5.4: Rankings of audible quality after reconstruction of shrimp clicks. Daubechies’ Least Asymmetric wavelets of order 8 and soft thresholding used

is 96% of the VisuShrink level whilst being 24% closer to the original wavelet decomposition.

The most important test is, of course, how good reconstructions from these thresholds sound in practice. Tables 5.3 and 5.4 show the rankings for sound reproduction quality as measured by the panel of musicians described earlier.

As can be seen, a threshold value of 0.1 was accepted as producing the best audio reproduction by the panel and this value is therefore used throughout this work. This result is not overly surprising, as VisuShrink reproductions were considered to be very poor in Chapter 4 and higher κ values produce results closer to VisuShrink.

It should be noted that whilst κ = 0.1 is appropriate for the particular classes of underwater signals analysed here, it is not necessarily the best choice for other types of signal. Indeed, it would be necessary to experiment in order to optimise the choice of threshold effect for any particular type of signal.

Significance level  0.1%  1%  5%  10%
Compression ratio 40.91% 40.39% 39.32% 38.83%
SS_{tot}(×10^{8})  27.55 26.46 22.97 21.09
%SS_{tot}  23.95% 23% 19.97% 18.33%
SS_{10}(×10^{8})  2.05  2.05  2.05  2.05
SS_{11}(×10^{8})  9.48  9.48  9.48  9.48
SS_{12}(×10^{8})  8.88  7.89  4.45  2.57
SS_{13}(×10^{8})  0.05  0.05  0.05  0.05
SS_{14}(×10^{8})  0.70  0.60  0.56  0.55
%SS_{10}  100% 100% 100% 100%
%SS_{11}  100% 100% 100% 100%
%SS_{12}  71.26% 63.33% 35.69% 20.65%
%SS_{13}  0.18% 0.18% 0.18% 0.18%
%SS_{14}  1.26% 1.08%  1%  0.99%

Table 5.5: Values of summary measures after reconstruction of dolphins whistles. Daubechies’ Least Asymmetric wavelets of order 8 and soft thresholding used

5.3.2 The Significance Level

One part of the runs test, as defined in Chapter 4, is the significance level used to decide if a block of coefficients contains a significant run of data. This significance level was arbitrarily set at 5% for the results given in Chapter 4 but is this the optimal choice?

Although a significance level of 5% for general statistical tests can be considered a reasonable default, there is an opportunity with the experiments undertaken in this work to discover if this is the best setting for underwater sound data of the type analysed here.

Accordingly, the summary measures used previously for various values of the significance level are shown in Tables 5.5 and 5.6.

These results were produced with a block size of 0.00625 seconds and a κ value of 0.1. This block size limits RunsThresh to levels above 11. Equally, as explained in Chapter 4, levels below 5 are not thresholded at all.

Again, the shrimp click data is much closer to the VisuShrink results than the dolphin whistle data. As expected, decreasing the significance level makes RunsThresh closer to VisuShrink. However, even with a significance

Significance level  0.1%  1%  5%  10%
Compression ratio 70.09% 67.36% 65.18% 63.38%
SS_{tot}(×10^{8})  16.04 15.25 14.76 14.32
%SS_{tot}  58.24% 55.37% 53.59% 51.98%
SS_{10}(×10^{8})  1.61  1.61  1.61  1.61
SS_{11}(×10^{8})  3.29  3.29  3.29  3.29
SS_{12}(×10^{8})  5.25  5.25  5.22  5.06
SS_{13}(×10^{8})  0.84  0.22  0.05  0.05
SS_{14}(×10^{8})  3.27  3.10  2.81  2.53
%SS_{10}  100% 100% 100% 100%
%SS_{11}  100% 100% 100% 100%
%SS_{12}  100% 100% 99.4% 96.38%
%SS_{13}  7.09% 1.84% 0.41% 0.41%
%SS_{14}  87.31% 82.81% 75.08% 67.49%

Table 5.6: Values of summary measures after reconstruction of shrimp clicks. Daubechies’ Least Asymmetric wavelets of order 8 and soft thresholding used

level as low as 0.1% there is still a considerable difference between the two threshold methods. There are not large differences in the compression ratio between the different significance levels which implies that it makes little difference which method is chosen, results which are borne out by the panel test below.

It should be noted that level 13 shows the most positive runs results, with the majority of the significance levels producing the same results. In general, there is very little difference between the results for significance levels of 1% and 5% with only level 12 of the dolphin data showing a major change. This implies that there should not be a major difference in reproduction quality between these significance levels.

The major test is, of course, how good the reconstructions sound. Accordingly, the experiments performed in the preceding section were repeated using this range of significance levels. Tables 5.7 and 5.8 show the results for dolphins and shrimps.

As can be seen, there was a general acceptance of 5% as an appropriate value, although 1% also performed reasonably. This is not too surprising as the other values made the sound far closer to the VisuShrink representation Significance

level (%)  Subject 1 Subject 2 Subject 3 Subject 4 Average
10%  3  2  3  3  3
5%  2  1  1  2  1
1%  1  3  2  1  2
0.1%  4  4  4  4  4

Table 5.7: Rankings of audible quality after reconstruction of dolphins whistles. Daubechies’ Least Asymmetric wavelets of order 8 and soft thresholding used

Significance

level (%)  Subject 1 Subject 2 Subject 3 Subject 4 Average
10%  3  3  2  3  3
5%  2  1  3  1  2
1%  1  2  1  2  1
0.1%  4  4  4  4  4

Table 5.8: Rankings of audible quality after reconstruction of shrimp clicks. Daubechies’ Least Asymmetric wavelets of order 8 and soft thresholding used

which was not felt to be very good by the panel in the results presented previously. The similarity between these two levels was also predicted by the summary measure results.

It should be made clear that in this test, the panel found it far more difficult to rank the sounds. Indeed, in the case of the shrimp clicks the average rankings show that 1% was the preferred choice. Generally, the panel felt that the different significance levels did not have a large impact on the sound quality, certainly the threshold level was held to be far more useful in producing different quality sounds.

As a result of these results, it was felt that there was no need to vary both significance level and κ simultaneously as the effect of changing κ would swamp any change in significance level.

5.3.3 The Block Size

The final value which could affect the results of RunsThreshing is the number of coefficients in each block which are tested. The value used in all the

Block size (s)  0.003125 0.00625 0.0125 0.025
Compression ratio 41.63% 39.32% 36.95% 36.22%
SS_{tot}(×10^{8}) 31.16 22.97 14.20 13.33
%SS_{tot}  27.09% 19.97% 12.34% 11.59%
SS_{10}(×10^{8})  2.05  2.05  2.05  0.17
SS_{11}(×10^{8})  9.48  9.48  3.77  3.32
SS_{12}(×10^{8})  12.46  4.45  1.39  1.31
SS_{13}(×10^{8})  0.08  0.05  0.05  0.05
SS_{14}(×10^{8})  0.71  0.56  0.55  0.55
%SS_{10}  100%  100% 100% 8.35%
%SS_{11}  100%  100% 39.76% 35.02%
%SS_{12}  100% 35.69% 11.13% 10.54%
%SS_{13}  0.27% 0.18% 0.18% 0.18%
%SS_{14}  1.27%  1%  0.99% 0.99%

Table 5.9: Values of summary measures after reconstruction of dolphins whistle. Daubechies’ Least Asymmetric wavelets of order 8 and soft thresholding used

experiments throughout this work is equivalent to a time of 0.00625s which is the smallest time period used by the Defence Research Agency. It therefore seems sensible that we follow the same conventions.

However, although the main set of results presented in Chapters 4 and 6 will use this window size for all such tests, we will examine here the effect of using differing block sizes. Results for the summary measures described earlier are shown in Tables 5.9 and 5.10.

These results were produced with a κ value of 0.1 and a significance level of 5%.

Shortening the block size brings the results closer to VisuShrink as fewer coefficients can be runs tested. However, the results are very similar at all lengths of block size. Again, this implies that the block size is not as significant as the κ level.

Again, level 13 shows the most positive runs results, with the majority of block sizes producing the same results. In order to determine the quality of reproduction, we again perform the panel test. The results are shown in Tables 5.11 and 5.12.

Block size (s)  0.003125 0.00625 0.0125 0.025
Compression ratio 68.17% 65.18% 63.28% 60.3%
SS_{tot}(×10^{8}) 15.74 14.76 14.49 13.69
%SS_{tot}  57.14% 53.59% 52.6% 49.69%
SS_{10}(×10^{8})  1.61  1.61  1.61  1.53
SS_{11}(×10^{8})  3.29  3.29  3.23  2.91
SS_{12}(×10^{8})  5.25  5.22  5.25  5.25
SS_{13}(×10^{8})  0.78  0.05  0.05  0.05
SS_{14}(×10^{8})  3.03  2.81  2.56  2.16
%SS_{10}  100%  100% 100%  95%
%SS_{11}  100%  100% 98.31% 88.45%
%SS_{12}  100%  99.4% 100% 100%
%SS_{13}  6.53% 0.41% 0.41% 0.41%
%SS_{14}  80.95% 75.08% 68.47% 57.81%

Table 5.10: Values of summary measures after reconstruction of shrimp clicks. Daubechies’ Least Asymmetric wavelets of order 8 and soft thresholding used

Block size (s) Subject 1 Subject 2 Subject 3 Subject 4 Average
0.003125  3  3  3  4  3
0.00625  1  2  2  1  1
0.0125  2  1  1  2  1
0.025  4  4  4  3  4

Table 5.11: Rankings of audible quality after reconstruction of dolphins whistles. Daubechies’ Least Asymmetric wavelets of order 8 and soft thresholding used

These results seem to indicate that the choice of block size between 0.00625 and 0.0125 seconds is difficult. This implies that, within a certain range, changing the block size does not produce a large effect. Outside this range, the block size will have an effect although it is unlikely to be positive. Too small a block size will obviously mean that fewer coefficients are RunsThreshed as fewer levels will contain enough values to be runs tested. Similarly, larger values beyond a certain point will not add anything as the levels with fewest coefficients are not thresholded in any event.

Either 0.00625 or 0.0125 would have been a good choice for the underwater data but as has been stated, external circumstances make us choose 0.00625 seconds which is 256 data points in length.

Block size (s) Subject 1 Subject 2 Subject 3 Subject 4 Average
0.003125  4  3  4  4  4
0.00625  1  2  1  3  2
0.0125  2  1  2  1  1
0.025  3  4  3  2  3

Table 5.12: Rankings of audible quality after reconstruction of shrimp clicks. Daubechies’ Least Asymmetric wavelets of order 8 and soft thresholding used

5.3.4 Conclusions

These results have shown that κ = 0.1, a significance level of 5% and a block size of 0.00625s are reasonable and sensible defaults with such data. Therefore, all the results throughout this work use these values.

This chapter has also shown that RunsThresh does behave significantly differently from VisuShrink and from VisuShrink with a globally reduced threshold. This implies that RunsThresh is a genuinely different and, as was shown in Chapter 4, useful threshold technique for the analysis of underwater sound data.

These results therefore confirm our conviction that using RunsThresh as a basis to build the remainder of our analysis upon is reasonable and sensible. The remainder of this work will therefore concentrate on using RunsThresh in preference to other threshold techniques discussed.

Signal Detection

6.1 Introduction

Whilst Chapter 4 was concerned with data compression and noise reduction, this chapter is concerned with identifying the actual positions of the signals, that is the points in the sound where the signals begin and end.

The wavelet transform allows one to examine the sound at different scales, much as the Fourier transform allows one to analyse data at different frequencies. As one of the unique identifiers for different signals is the scale (analogous to frequency) at which certain changes take place, it is perfectly possible to use the wavelet transform for signal identification.

For the analysis in this work however, the Fourier decomposition of a sound has one advantage over the wavelet decomposition. In a Fourier representation, each frequency band has exactly the same number of observations as every other frequency band, allowing direct comparisons across frequencies by way of standard statistical techniques (such as principal component analysis, projection pursuit or even Chernoff’s faces for a graphical representation). The wavelet representation however, due to the pyramid algorithm described in Chapter 2, doubles the number of coefficients in each successive scale band. This obviously makes simple comparisons more difficult when using the standard wavelet decomposition.

As we have used the wavelet transform for data compression and noise reduction, it will be highly advantageous to use it for signal identification as well. For example, the ability to perform only one transform for all three processes will save time. In order to do this some adjustment to the wavelet transform is needed in order to obtain a rectangular coefficient matrix rather than a triangular one.

There are two obvious ways to approach this problem: one is to alter the transform itself in some way so that it produces equal numbers of coefficients at each scale; the other is to post-process the coefficients at each scale so that there are equal numbers of summary measures at each level.

Here, the second of these methods is the one adopted and this is described in some detail below. This approach was chosen as it will also produce a reasonable reduction in the number of coefficients to be analysed which, given the amount of data, is an important side effect. Nason and Silverman (1995) describe a possible way of implementing the first approach via the Stationary Wavelet Transform, developed by Pesquest et al. (1994).

Once the wavelet decomposition has been chosen as the starting point for the signal identification, there is still a wide choice of techniques available to actually identify parts of the decomposition that represent signals.

For this work, the choice is restricted somewhat however by the fact that there is no possibility of producing “training sets” of known signals for comparison purposes. This excludes such widely known methods as neural networks — Ripley (1994) provides an introduction — and Kalman filtering, which is a commonly used signal detection method in engineering — McCausland (1969) gives a good introduction to this field.

It was therefore decided that, as the method should be non-parametric and indeed data driven, some form of multivariate density estimation might be the best choice for actually describing the noise in the data.

Thus, this chapter begins by describing the method of post-processing the wavelet decomposition to produce equal numbers of coefficients at each scale level, and follows this with a discussion of the different kernel smoothing techniques that could be applied, including the chosen method. Next, the technique used to actually identify signals given these methods is shown and finally the results of this method of analysis for the underwater sound data are presented, comparing these with a Fourier based technique.

6.2 Squaring the Pyramid

The first task in the identification sequence, as shown in Figure 6.1 is to transform the wavelet decomposition pyramid into a standard matrix. Although the method chosen here is based on coefficient summarisation, a brief note on the alternative approach using the Stationary Wavelet Transform (SWT) is first given.

Figure 6.1: The process of sound identification.

6.2.1 The Stationary Wavelet Transform

The Stationary Wavelet Transform was developed by Pesquest et al. (1994)

— from an engineering perspective — and described in a statistical context by Nason and Silverman (1995).

This modification attempts to square the pyramid by filling in the missing coefficients in lower levels of the wavelet decomposition. This obviously introduces some redundancy into the results, as the transform contains N/2 coefficients at each level, but the resulting diagrams are easier to analyse by eye. Figure 6.2 shows the normal DWT and the SWT for a standard Doppler signal (as described in Chapter 4).

The detailed mechanics of the Stationary Wavelet Transform will not be given here but briefly, the SWT comes in two forms: a time-ordered version useful for time series analysis and a packet-ordered version useful for curve estimation. Obviously, the time-ordered version is the most useful for this work.

(a) Discrete (b) Stationary

Figure 6.2: Two wavelet decompositions of a Doppler signal.

The time-ordered SWT simply collects all possible (dyadic) shifts of the discrete wavelet transform into one tableau. Unlike the original DWT, this means that the “SWT coefficients of translated data are identical to the translated coefficients of the SWT of the original data” (Nason 1997). This transform is O(log N) slower than the standard DWT and is achieved, in principle, by not performing the decimation step in the standard transform.

It is worth considering here how this method could be used in signal identification. Figure 6.2 shows that with each level being filled in, the coefficients in any one level lend themselves to density estimation using a univariate kernel. This would suggest a series of multiple estimates (one for each level) with any deviation from the density being flagged as a signal. A threshold could then be introduced which would trigger identification if a certain number of levels showed a signal.

Such an approach is obviously worth investigating and should prove useful for certain applications, indeed Nason and Silverman (1995) use a kernel smoothing approach to analyse shock impact data. However, as the approach drastically increases the number of coefficients from each wavelet decomposition, and as the underwater sound data sets are already massive, this work follows an alternative, summarisation based, approach which will now be examined in detail.

6.2.2 Summarising the DWT

The method developed here is a natural extension to the ideas presented in Chapter 4 with RunsThresh. There, it was required to partition the decomposition into coefficient blocks in order to perform a runs-test on them for the thresholding and similarly, blocks of coefficients are needed in order to summarise the coefficients.

The major difference between this approach and the SWT is as follows. Whilst the SWT works for all levels of coefficients, the amount of data generated is actually vastly higher than with the plain DWT. A DWT of 256 data points will produce 256 coefficients, whilst a SWT will produce 1024 coefficients. This obviously means that more computing power is needed “downstream” to analyse the matrix. By contrast, summarising the DWT coefficients produces an element of data reduction, easing the processing burden for the actual analysis, but means that not all scale levels can be included, as lower frequency levels have too few coefficients to be in any of the summarised blocks. This was not a problem with thresholding, as the lower frequencies tend to be excluded from thresholding anyway.

To a large extent, the choice of method depends on the intended final stage of the identification problem. As this work is going to concentrate on using kernel based methods for the signal identification, and as these methods are highly computer-intensive, it makes sense to use the summarisation technique as this will reduce the amount of work that needs to be performed, thus speeding up the whole detection procedure. For other applications, or other final identification methods which need fewer resources, the SWT may be the preferred choice.

Choice of Summary Measures

Given that the coefficients in each level are going to be summarised, there are two pertinent questions that need to be answered — what size should the summary windows be and how are the coefficients to be summarised?

The first of these questions is the easier to answer as there is a natural limit to the window size imposed by the choice of scale levels that are to be summarised. As the frequencies which are required decrease, a longer window is needed so that these low-scale coefficients can be included. There is a secondary factor, which is the length of signals that are required to be detected. As signal identification is to be based on the summarised coefficients, it would be unwise to have windows that are much longer than the anticipated signals as this will avoid having multiple sounds in the same window which would then not be detectable as individual events.

The second question is more difficult to answer, because there are far more options. Fristrup and Walkins (1993) have enumerated a number of possible variables that change over the course of a sound and could be used to distinguish signals from background noise. Obviously, with such a vast number of potential choices, it is best to start with the simplest possibilities before trying more complex, and time consuming, ideas.

The simplest possible measure to describe a sound is the energy. The total energy of a discrete sound is the sum of the squares of the data values. As the DWT preserves the total energy, it follows that the sum of the squares of the wavelet coefficients is also equal to the total energy, and in fact in each scale level the energy is preserved so the mean sum of squared coefficients in each level is the mean energy at that scale.

With this in mind, the first chosen summary measure was the mean square energy for windows at each frequency level. The preservation of energy by the transform means that although lower scale levels have fewer coefficients in such windows, the values of these coefficients are higher so that the mean is the same at each scale level — for a sound that has equal energy across the whole frequency range of course.

A formal definition of a new coefficient, e_{i}, at level j is therefore:

e_{ij} = \frac{1}{2^{j}} ∑_{k=1}^{I/2^{ln(n)−j}}(w_{j,k})^{2}

where k runs over the coefficients in the window and I is the number of windows selected.

These coefficients form a matrix, e_{i,j}, over the levels j which are less than the complete range of scale levels available from the wavelet decomposition. Now, while this is an unfortunate side-effect of this method of “squaring the pyramid”, it is not a problem with the data sets examined in this work. Underwater sounds are, as has been shown in Chapter 3, broad band, with a typical range of 0–16 kHz showing some activity. However, initial examination of the data shows that differences in the signals of interest are most marked in the higher frequency bands, therefore the loss of the lower frequency elements of the wavelet decomposition is not the problem that it may have been for other types of data.

For the particular sounds used in this work, it was decided that the highest eight scale levels exhibited the most important multivariate behaviour for differentiating between the different signals. As the shortest signal that detection is required for is a click, as shown in Chapter 3, and this is measured in hundredths of a second, it was further decided that a window of 256 data points, or 0.00625s, would be a suitable size, ensuring that the chance of two clicks occurring in the same summary window is minimised.

6.3 Signal Identification

Preliminary investigations, as shown in Bailey et al. (1996), suggested that sound identification could be achieved via detection of changes in the joint behaviour of the mean square energy at several scale levels. To determine if the joint behaviour is changing, a reference model of the joint distribution during periods of noise is required along with a method of detecting multivariate outliers from this distribution.

If some simple assumptions about the multivariate distribution during periods of noise could be made, such as multivariate normality, then the parameters of this distribution could be readily estimated and various techniques would be available to detect outliers. Penny (1996), for example, uses the Mahalanobis distance to achieve this.

The complex and, more importantly, changing nature of the background underwater environment precludes any such simple assumptions however. The use of a non-parametric and flexible data dependent approach to the estimation of the multivariate probability density would, in fact, be superior. Kernel probability density estimation was selected as a suitable nonparametric method for this method.

After forming the summary matrix e_{i,j}, described above, the multivariate behaviour across levels may then be subjected to analysis. This may be contrasted with the approach taken by Nason and Silverman (1995) where they use a kernel approach in each level but where they do not look at behaviour across levels.

The kernel smoothing technique for probability density estimation will now be explained in detail. After this, a description of how this estimate is used to identify signals will be given.

6.3.1 Kernel Density Estimation

Several texts have details of the basic density estimation approaches possible, including Silverman (1986), Wand and Jones (1995) and Scott (1992), all of which give a good general introduction to this area of statistics.

Silverman (1986) gives a clear description of the kernel method in several dimensions, and much of the following description is based on remarks in this text.

Kernel density estimation aims to form an estimate of the probability density function of a supplied data set. In this case the data set is multidimensional and hence a multivariate probability density will be found.

Now,as we have defined e_{i} as a series of points in j dimensional space, corresponding to the multivariate summary measures described in the preceding section, then the kernel estimator can be thought of as a series of ‘bumps’ centred at each e_{i}. Therefore, the multivariate kernel density estimator with kernel K and window width h for n data vectors is defined by

f̂(e) = \frac{1}{nh^{j}} ∑_{i=1}^{n} K\{\frac{1}{h}(e − e_{i})\}

The two important factors in this definition are the kernel function and the bandwidth (or window width). Each of these will now be discussed.

The Kernel Function

The j dimensional kernel function, K(e), is a weighting function, satisfying

∫_{ℜ^{j}} K(e)de = 1.

Although K could be any shape, chosen from the set of all functions that meet this condition, a spherically symmetrical unimodal probability density function is generally used. A common example is the standard multivariate normal density, given by

K(e) = (2π)^{−j/2} exp(−\frac{e′e}{2}).

Various other functions such as the multivariate Epanechnikov kernel and, for the two dimensional case, K(e) = 3π^{−1}(1−e′e)^{2} are also used, and indeed, allow faster computation than the standard multivariate normal.

As the data are not necessarily spherical, it would be wise to pre-scale the data, as is done in most multivariate statistical approaches. Although this could be performed prior to the density estimation, Fukunaga (1972) has developed an attractive method which combines scaling and kernel density estimation in one process.

The idea, shown in Figure 6.3, is to use a kernel that is the same shape as the data itself. Linearly transforming the data to have a unit covariance matrix then smoothing using a radially symmetrical kernel before finally transforming back is equivalent to using such a kernel.

Figure 6.3: Kernel and data shapes if the Fukunaga method is used, after Fukunaga (1972).

This process is actually equivalent to simply using a density estimate of the form

f̂(e) = \frac{(det S)^{− \frac{1}{2}}}{nh^{j}} ∑_{i=1}^{n} k\{h^{−2}(e − e_{i})′S^{−1}(e − e_{i})\}

where k is given by

k(e′e) = K(e)

and S is the sample covariance matrix, with K being the normal kernel, k(u) = (2π)^{−j/2} exp(−\frac{1}{2}u). It is may be advisable for S to be robustly estimated, using, for example, the method suggested in Tukey and Tukey (1981) section 10.2.2, especially if there is the possibility of the data containing extreme outliers.

The summary measures used in this analysis are highly likely to be on vastly different scales, necessitating the use of scaling prior to analysis. This restriction makes the Fukunaga kernel a sensible choice, allowing the data to be used in the density estimation without any prior processing.

As noted in Wand and Jones (1995) however, there is no theoretical justification for sphering data that is not multivariate normal. The data used here are highly unlikely to be multivariate normal and indeed the density function may be changing over time. Analysis of our data sets show that scaling the data is appropriate, even if only as a computational aid. This requirement for scaling the summary measure data, together with the results obtained, have strongly influenced the decision to use the Fukunaga kernel method in the remainder of this work.

The probability density estimate used for the remainder of this work is therefore of the form:

f̂(e) = \frac{(det S)^{− \frac{1}{2}}}{(2π)^{\frac{j}{2}} nh^{j}} ∑_{i=1}^{n} exp\{−\frac{1}{2h^{2}} (e − e_{i})′S^{−1}(e − e_{i})\}

It should also be noted that Friedman et al. (1975) have developed a computationally faster procedure for evaluating the kernel at irregular points (as is the case is here). This involves first ordering the data along the coordinate in which the spread is greatest and calculating the kernel ignoring observations which cannot contribute to the value at the point being looked at. This can obviously decrease the computation time significantly.

The Bandwidth

The other important factor in the density estimate is the bandwidth, h. Although this is shown as a single variable, it could be replaced with a matrix, h, with different bandwidths for different variables. In practice this is found not to be necessary as using a single bandwidth with scaled data gives similar results to the use of a bandwidth matrix and unscaled variables. Such scaling is performed by the Fukunaga method and is, of course, the main reason for performing the density estimation in this manner.

The question of how the bandwidth is chosen is also important. For a single bandwidth or multiple bandwidths there are two main choices, adaptive (data-driven) methods or “optimal” fixed methods. In this discussion, the optimal choice of a single bandwidth is considered. These ideas can be extended to the case of multiple bandwidths however.

The question is therefore how to find this optimum bandwidth, h_{opt}. As mentioned, there are two general approaches that can be used, fixed value type methods and data-driven methods. In the former, the bandwidth is selected according to some pre-specified rules or heuristics whereas in the latter, h_{opt} is determined by using some function of the data. The main differences are speed and reliability, with data-driven approaches being slower but more accurate, especially for data sets that depart from multivariate normality.

As an example of a fixed value method, consider the case when f has bounded and continuous second derivatives. It can be shown that the optimal bandwidth converges to zero as the number of data points increase, but at an extremely slow rate. It can also be shown that the bandwidth depends on the unknown density. Now, Silverman (1986) suggests that one should initially find the optimal value for h when f is a standard density. In the case that a standard j-variate normal density is assumed, he shows that h_{opt} is given by

h_{opt} = A(K)n^{− \frac{1}{j+4}}

where the constant A(K) varies with the choice of kernel. For the standard multivariate normal kernel used in this work, A(K) is {4/(2j + 1)}^{1/(j+4)} (which is 0.96 for the two dimensional case).

As the Fukunaga method is being used, this value is directly usable as the optimal bandwidth, for non-scaled data, the data should be scaled by a factor σ and the optimal bandwidth should be set to σh_{opt} where σ^{2} is, for example, the average marginal variance.

As is noted in the above text, this method “does give a quick and easy choice of at least an initial value for the window width” but may not be optimal. It is however, obviously quick to calculate, something that is not the case for the best data-driven methods.

So which data-driven techniques would be suitable to produce estimates of the optimal bandwidth? Both Wand and Jones (1995) and Silverman (1986) suggest using some form of cross-validation to discover the optimal bandwidth for any given set of data. Although Silverman makes the case for using a least-squares cross-validation method, he acknowledges that the likelihood method, as advocated by, for example, Habbema et al. (1974), produces satisfying results if using an unbounded kernel (as is the case in this work) and for data with no extreme outliers — in the case considered here of course, the data is, by definition, free of outliers.

This work implements the likelihood cross-validation method for choosing an optimal bandwidth, h_{opt} when one has a matrix of n multivariate observations, x_{1}, …, x_{n} from which one is estimating a density.

The reasoning behind this approach is as follows, after Bowman (1984). For the entire n point data set used to form the kernel density estimate, an observation has obviously occurred at every point. Hence, the likelihood of the entire set is approximately

∏_{t=1}^{n} f̂(x_{t})

where the x_{i} are fixed and the total product is a variable function of the bandwidth, h. An obvious approach to choosing the optimal h is therefore to choose a bandwidth which maximises this expression.

This simple approach has one problem however, for such an expression there is always a global maximum at h = 0. An adaption is required to circumvent this maximum. If the approximate likelihood is expressed as

∏_{t=1}^{n} f̂_{t}(x_{t})

where f̂_{t}(x_{t}) refers to the kernel estimate using all observations excluding x_{t}, the log-likelihood is obviously

∑_{t=1}^{n} log(f̂_{t}(x_{t}))

in this case. This is the likelihood cross-validation function.

This function, suggested by Duin (1976) and Habbema et al. (1974), does rely on various assumptions, some of which are overcome if K has unbounded support —fortunately, that is the case here. Bowman (1984) has suggested a slight modification that is less sensitive to large numbers of extreme outliers.

Although this presents a method to determine the optimum bandwidth, a search still needs to be performed to enable the actual value of h_{opt} to be discerned. The search method chosen in this case was a golden ratio search, as explained in Press et al. (1992). This is a reliable method of finding a local optimum of, in this case, the bandwidth.

To perform the search, two fixed bandwidths are used and a third is then found. The fixed points are 0, where the score is set to −∞ to force the search to exclude a bandwidth of 0, and the previously found optimum bandwidth, h_{opt}, with its corresponding score. For the first run of the routine, the value of h_{opt} is set to

h_{opt} = (\frac{4}{n(2j + 1)})^{\frac{1}{j+4}}

as explained above.

This gives a lower bound and a mid-point, the upper bound for the search is then found by trying the two points that force the mid-point to split the range into a golden section, that is where the sections are split in the proportions 0.61803399 and 0.3819661 (these are the golden ratios). Obviously, the point that is tried first is the one which would be closest to the mid-point, to keep the range as narrow as possible. If both of these bandwidths fail to form a boundary — the score function does not turn — then the upper bound is continuously doubled until a limiting bound is found.

The standard golden section search is then performed, with new bandwidths being tried in such a way as to retain the golden ratio split between the sections. Eventually, this will iterate down to a single bandwidth which is then declared to be the optimum and used in the calculation of the kernel density estimate.

This process, although a robust method for finding the optimal bandwidth, is slow for two reasons. Firstly, it involves multiple calculations of the kernel values at each stage, although this would be a common feature of any search algorithm. More importantly, it is inefficient, that is, the search does not converge very fast. It is however, a necessary process. The optimal bandwidth can be wildly different from the value A(K)n^{−1/(j+4)} suggested by Silverman (1986), although this is not surprising, given the conditions under which this value is valid.

It should of course be noted that this is just one of many data-driven methods for optimal bandwidth selection. Other techniques include least squares cross-validation, various forms of “plug-in” rules, smoothed crossvalidation and biased cross-validation. Further works that discuss this issue include Wand and Jones (1994), Park and Marron (1990) and Park and Turlach (1992), whilst Wand and Jones (1995) give a comprehensive overview.

6.3.2 Identifying Sounds Using Probability Density Estimates

Following through the signal identification process so far (from Figure 6.1), the wavelet decomposition has been summarised, in terms of the energy in each scale level, and the methods discussed in the previous section have enabled a joint probability density to be estimated.

It should be noted at this point that it may not be necessary, or even desirable, to form the density estimate based on all the summary variables.

Investigation has shown that a subset of these variables are sufficient to describe many of our sounds in such a way that unique events are differentiated by the subset of variables. In such a case it would be desirable to use the smallest subset possible for computational reasons, the speed of optimum bandwidth and probability density estimation both being functions of the dimensionality of the data being analysed. This topic is examined further below, after an explanation of the rest of the identification process.

Given the impossibility of producing training sets, due to the restriction that the method must be able to cope with signals that have never been previously encountered, it is essential to impose one restriction on the application of the proposed method. It is imperative, and indeed a necessary condition, that any sound to be analysed using the following procedure must have a section at the start that is known to be background noise, contaminated only by events that are not of interest. The required length of this section of noise depends on several factors which will be described below.

Given a matrix of summary coefficients, e, with an initial number of observations representing noise, the probability density can then be estimated for these initial observations. This then gives a statistical description of the noise against which new observations can be compared. The general procedure is to take each new observation in order and compare this with the density estimate. If the observation is an outlier from the density then it is flagged as a signal of interest, otherwise it is declared to be noise.

If the observation is noise then the oldest observation involved in the estimate of the density is removed and the new observation is included. The optimal bandwidth for the new data set is then found and the process repeated. There are various methods for deciding if the observation is an outlier from the multivariate density estimate and these are described below.

An optional extra smoothing effect can be introduced by using a buffer to lag the updating of the density estimate. If an observation is found to be noise then, instead of being added to the data that is kernel smoothed, it is added to a buffer of a length suitable for the data. If the buffer is full, then the oldest buffer value is added to the set of observations used to estimate the probability density and removed from the buffer. If the buffer is not full, the density estimate is not updated. In effect, this enforces a delay before the density estimate is updated, which can help ensure that signals are found if they gradually rise from the background noise rather than having a more abrupt start. It also helps to minimise a false negative effect being spread forwards through time —a false negative will be discounted by the length of the buffer, minimising its effects.

It should be noted that the estimate of the background noise density adapts recursively over time as the data on which the it is built constantly updates. This allows any gradual changes in the background underwater environment that occur to be represented as they happen.

Some of the factors that were mentioned above will now be examined in more detail.

Number of Variables to Analyse

As was noted above, there is quite a convincing case to be made for using the minimum subset of variables necessary in the analysis. This case is based on computational speed and the fact that smaller dimensionality requires substantially fewer initial noise values for accurate estimates of density values (as shown below). Given that this is the case, a decision needs to be made as to how large a subset should be examined and what variables this subset will consist of.

Figure 6.4 shows a plot of the eight summarised variables (being wavelet decomposition levels 7–14) used in this work for a section of rain falling, a dolphin whistle, a seal bark and shrimp clicks.

As can be seen, the differing signals all exhibit different characteristics. However, the graphs are very confused and this difference of behaviour can be seen by examining a subset of these levels. In this case, if levels 11–13 are examined (see Figure 6.5), the same differences are present and can be seen much more clearly.

These variables are created from an original wavelet decomposition of 32768 data points (0.8s) and each window is 256 data points in length (0.00625s). These variables are therefore equivalent to scale bands equat-

(a) Rain  (b) Dolphin Whistle

(c) Seal Bark  (d) Shrimp Clicks

Figure 6.4: The total energy by variable for various signals.

ing to frequencies of 2.5, 5 and 10 kHz.

It should also be noted that the behaviour is only different for each type of sound if the variables representing all three scale levels are jointly considered — the multivariate behaviour is the key to identifying each type of sound rather than simply one dimensional behaviour, hence the necessity of a multivariate kernel estimator throughout this work.

Although the underwater sound recorded in the conditions and by the method specified here exhibits distinctive behaviour in this set of frequency ranges, that is obviously no guarantee that these variables will be a good

(c) Seal Bark  (d) Shrimp Clicks

Figure 6.5: The total energy by a subset of variables for various signals.

choice for every data set; indeed different recording circumstances even for this data may require a different subset of variables.

There is no easy method of deciding upon a subset for unknown data. If computational ability is not a problem then the complete set of variables can always be analysed or, if some a-priori information about the data set is known, as is the case here, then a subset can be selected by this method. Indeed, analysing the wavelet decompositions to gain some knowledge about the signals being sought allowed a subset of variables to be selected with some speed for the underwater data.

Finally, a method involving cross-validation techniques could be used with a suitable set of known data to achieve some form of automatic selection of scale levels but that is a topic for future research which will be discussed in chapter 7.

Length of Initial Noise

In the case of univariate space, it is reasonably easy to calculate what sample size is required to estimate accurately the probability density at any one point. In the multivariate case however, the problem is more difficult. Silverman (1986) has produced a table (Table 4.2, p94) of required sample size against number of dimensions for accurate estimation of the density at a point in the main area of the probability density, as opposed to in the tail. Intuitively, for accurate estimation of tail values, a larger sample size would be expected to be required. As this process relies on detecting outliers this is quite an important point.

Dimensionality Required Sample Size Time (s)
1  4  0.025
2  19  0.120
3  67  0.420
4  223  1.40
5  768  4.80
6  2790  17.0
8  43700  270
10  842000  5300

Table 6.1: Sample size required to produce acceptable mean square errors in the estimation of the density at 0, after Silverman (1986).

Table 6.1 gives an indication of just how quickly the required size of the initial noise section increases with dimensionality. Assuming windows with a length 0.00625s are used for the underwater sound the value for 10 dimensions corresponds to a time period of approximately an hour and a half.

As can be seen, required sample size increases rapidly with the number of dimensions. Furthermore, these figures are given for the point 0. As Silverman (1986) states, this means that the sizes given are likely to be an underestimate of the sample size required for accurate measurement in the tails of the distribution.

Obviously, this implies that the number of dimensions needs to be kept as low as possible, especially if the length of sound to be analysed is quite short.

In fact, if the sound to be analysed has been pre-recorded, rather than being analysed as it is captured, the nature of the recording could impose a limit on the number of dimensions due to the length of the initial noise. For instance, in the 12s of sound analysed later in this chapter, the initial noise is no more than 1s (or about 150 observations) of the total recording, meaning that at most three variables can be employed.

Methods of Outlier Detection

As will be recalled from the earlier discussion, the identification of a signal depends upon it being flagged as an outlier from the estimated probability density.

There is a variety of possible ways that can be used to decide if a point is an outlier from a multivariate distribution. Six methods of various complexity, two of which require a kernel value and five based on Mahalanobis distances, are considered here.

Consider a new observation x_{n+1} that is to be classified as noise or signal based upon the density estimate of the noise built from n immediately preceding points, x_{i}, i = 1, …, n. One of these five tests of whether x_{n+1} is an outlier from this current kernel density estimate depends solely on the value of the density estimate at x_{n+1}, four are based solely on the Mahalanobis distance and one combines both these measures.

The six outlier test statistics, with δ being the squared Mahalanobis distance, are:

• The kernel value at the new point, f̂(x_{n+1}).

• The squared Mahalanobis distance from the new point to the mean of the points in the density, δ(x_{n+1}, \overline{x}).

• The squared Mahalanobis distance from the new point to the mean of the points in the density, weighted by the kernel values at those points, δ(x_{n+1}, \overline{x_{w}}).

• The mean squared Mahalanobis distance from the new point to each of the points in the density estimate, \overline{δ(x_{n+1}, x_{t})}.

• The squared Mahalanobis distance, weighted by the kernel values, \frac{∑ f̂(x_{t})δ(x_{n+1},x_{t})}{∑ f̂(x_{t})}.

• The squared Mahalanobis distance from the new point to the point in the density with the largest value of f̂ (x_{m}), δ(x_{n+1}, x_{m}).

These test statistics are compared with similar statistics built upon the points in the density to decide if the points represented are, in fact, outliers. All these comparisons are adaptive, that is, they change over time with the underwater sound background — they are, in a sense, analogous to the crossvalidation method of bandwidth selection.

For the kernel value test statistic, the new point is declared to be an outlier if

f̂(x_{n+1}) < min (f̂_{1}(x_{1}), …, f̂_{n}(x_{n}))

The comparison for the second test statistic is given by

δ(x_{n+1}, \overline{x}) > max (δ(x_{1}, \overline{x}_{t≠1}), …, δ(x_{n}, \overline{x}_{t≠n}))

whilst for the third statistic, the comparison is given by

δ(x_{n+1}, \overline{x_{w}}) > max (δ(x_{1}, \overline{x_{w}}_{t≠1}), …, δ(x_{n}, \overline{x_{w}}_{t≠n}))

The mean squared Mahalanobis to the points in the density declares the new point to be an outlier if

\overline{δ(x_{n+1}, x)} > max (\overline{δ(x_{1}, x_{t≠1})}, …, \overline{δ(x_{n}, x_{t≠n})})

the weighted squared Mahalanobis distance test declares x_{n+1} to be an outlier if

\frac{∑_{i} f̂(x_{t})δ(x_{n+1}, x_{t})}{∑_{t} f̂(x_{t})} > max (\frac{∑_{t≠1} f̂_{1}(x_{t})δ(x_{1}, x_{t})}{∑_{t≠1} f̂_{1}(x_{t})}, …, \frac{∑_{t≠n} f̂_{n}(x_{t})δ(x_{n}, x_{t})}{∑_{t≠n} f̂_{n}(x_{t})})

and finally, for the last test statistic, the comparison is given by

δ(x_{n+1}, x_{m}) > max (δ(x_{1}, x_{m≠1}), …, δ(x_{n}, x_{m≠n}))

In all these comparisons, f̂(·) refers to the kernel density estimate based on x_{1}, …, x_{n}; f̂_{t}(·) refers to that based on x_{1}, …, x_{n} but excluding x_{t}; and where δ(·) again refers to the squared Mahalanobis distance between two points. It should be noted that f̂_{t}(·), t = 1, …, n, are automatically available from the last iteration of the most recent cross-validation used to optimise the kernel bandwidth.

The first of these tests represents a simple comparison of f̂(x_{n+1}) with n examples of effectively the same quantity for points known to be noise, these examples being obtained by cross-validation on x_{1}, …, x_{n}; this type of criterion is common in kernel classification methods as shown in, for example, Remme et al. (1980).

The penultimate test is possibly the most unusual. Whilst it is also based on the Mahalanobis distance, the measure compared in this case is the weighted average squared Mahalanobis distance of x_{n+1} from x_{i}, i = 1, …, n where the weights reflect the estimated probability of x_{i} occurring under the current noise model. It was felt that involving the estimated current noise model in this way would improve the probability of detecting outliers from an adaptive background.

The remaining tests are based on a similar cross-validation idea, but the measure compared in these cases are various adjustments of the standard Mahalanobis distance. These ideas are modifications of the standard Mahalanobis distance criterion commonly used for detecting single multivariate outliers, as used in some of the other tests described by, amongst others, Atkinson and Mulira (1993) and Penny (1996).

The choice of outlier detection method actually used depends to some extent on the data under examination, although pointers to good general methods are given in the results below. Comparisons of type I and II errors for all the methods described above are included, based on the “operator” data with known signal described in Chapter 3.

It is also conceivable that several tests could be run in parallel and some combination of the results used to decide if a new observation is sound or noise.

6.3.3 The Role of Thresholding in Identification

As was shown in Chapter 4, thresholding the data is a useful method of cleaning the sound by reducing (or eliminating) noise. This leads to the question of what role, if any, thresholding should play in signal identification.

Obviously, even if a threshold is not applied before signal identification takes place, it can still be applied afterwards to enhance the audible quality and reduce the storage requirements of the sound. This reduces the question to the simple one of whether the sound should be thresholded before the signal identification is attempted.

In theory, the effect of applying the threshold before any identification will be to reduce the general background noise towards zero, whilst leaving the signals unaffected. This should make the task of the identifier easier as the signals will be more distant outliers from the noise than they would otherwise be.

However, no threshold can be expected to differentiate perfectly between noise and signal and indeed, there is a possibility that parts of the signal would be thresholded out, as shown in Chapter 4. It is therefore likely that either not all of the noise will be reduced, or that some of the signal would be reduced or possibly even that both would occur.

The first case would result in some of the noise being more likely to be mis-classified as signal, increasing the type I error rate, whilst the second case would result in signals being more likely to be classified as noise, increasing the type II error rate. The third case would obviously result in an increase in both type I and type II errors.

Overall then, although there are theoretical benefits to be gained from pre-thresholding the data, practical constraints on the quality of thresholding possible mean that it is more wise only to allow thresholding after the signal identification has taken place.

6.4 Results and Comparisons

In order to examine the usefulness of the methods detailed above, they were all used to analyse two types of data, described in Chapter 3. The first was the “operator” data with known signal positions and the second was various underwater sounds with unknown signal positions.

These recordings were subjected to the procedure described above with a variety of choices for such options as the initial noise size, using the program KP, or Kernel Predictor, described in detail in appendix B.

As a useful comparison, the Defence Research Agency supplied a program, Spectool, to subject the recordings to a Fourier based analysis. This is more fully described below.

We first give a description of the Fourier approach, then present the results for the two methods, and finally give some comparisons.

6.4.1 The Fourier Technique

The Spectool program implements two Fourier segmentation based methods of detection. Both work on the spectrogram described in Chapter 2 and are variants of spectral detection methods.

In common with the kernel method described above, a section of known noise is required, so that some basic statistics can be calculated. Either of the two methods can then be used — these are energy and spectral detection. Such methods are standard approaches to the analysis of sound using Fourier transform based techniques, similar work has been done by Boashash and O’Shea (1990) for example.

Energy Detection

This method of sound identification calculates the total energy within a time window and compares this with its expected value. The window is identified as sound if the difference is above a preset threshold. That is, a detection is made in a window if

\frac{power − mean(power)}{σ} > α

where σ is the standard deviation of the power, with this and the mean being calculated for an initial section of sound known to be noise. The threshold value is given by α.

The choice of threshold is obviously an important factor in the accuracy of the detection results. As α is lowered both the correct detection results and the false alarm rates rise. Therefore, the threshold must be set at a high enough rate to minimise the errors whilst still detecting the sounds.

In practice, an intuitive approach is to set the threshold at the lowest level possible that gives no false positives in the initial section of sound known to be noise. This is not ideal, but given that there is no chance to use a crossvalidation approach within Spectool, it is a sensible choice — the noise section is, after all, known to be noise.

Spectral Detection

Whilst the energy detection assumes that the sound has broadly flat spectral characteristics, and is therefore a rather crude method, spectral detection looks for deviations from the norm at individual frequencies. This is more similar to, although as will be seen, less advanced than, the kernel method described previously.

The FFT is used to calculate the energy at a particular frequency, x, within a time window, t, and this is compared with the mean for that frequency bin. Two thresholds, α and β, are needed. The former is a frequency bin threshold whilst the latter is for the total number of frequency bins used, according to the following formula:

f(x) = \begin{cases}1 & if \frac{power(x_{t})−mean(power(x))}{σ_{x}} > α \\ 0 & otherwise \end{cases}

where σ_{x} and mean(power(x)) are the standard deviation and mean of the power at frequency x, calculated for a section of sound known to be noise. Then

s = ∑_{x=0}^{n} f(x)

n being the number of frequency bins. A detection is flagged if

s > β.

Obviously, β needs to be reasonably high to prevent too many false alarms, but this will make the method biased towards broad band signals such as clicks, or signals with distinct harmonics, such as whistles. Given a suitable range for β, 15–25 is typical, then α can be chosen in the same manner as for the energy detection, being set as low as possible such that the noise section has no false positives found.

Results

Before presenting the results, a few words are necessary to explain the requirements for this work. The Defence Research Agency has confirmed the author’s belief that failing to detect signals is a far more serious problem than that of classifying noise as signal. This is so much the case that the DRA have specified that a false positive (or type I error) rate of anything up to 20% is entirely reasonable as long as this helps to reduce the false negative (or type II error) rate.

The analysis of error rates on actual underwater sound is extremely difficult as the true locations of signals are unknown to any degree of accuracy. As explained in Chapter 3, this has led to the use of an artificial data set, containing the word “operator” superimposed over “pink” noise with similar characteristics to the underwater environment. This data set has been created at a range of signal to noise ratios (SNRs) to allow results to be obtained in a variety of conditions. With an SNR of 0 dB the noise and signal have equal power and at -10 dB the signal is virtually indistinguishable by the human ear.

The initial investigation examined the choice of outlier detection measure. Tables 6.2 through 6.7 show the error rates for different SNRs and outlier detection measures. In all these cases, the density estimate was built on the initial 30 data points of the sound, with three summary variables being used for the analysis, corresponding to levels 11–13 of the wavelet decomposition.

SNR Type I (%) Type II (%) Sound identified (%)
+5  4.7  0.0  100.0
+1  11.8  0.0  100.0
0  2.4  7.7  92.3
-1  2.4  7.7  92.3
-5  7.1  38.5  61.5
-8  1.2  61.5  38.5
-10  9.4  46.2  53.8

Table 6.2: Error rates for the operator data using f̂(x_{n+1})

SNR Type I (%) Type II (%) Sound identified (%)
+5  4.7  0.0  100.0
+1  11.8  7.7  92.3
0  2.4  7.7  92.3
-1  2.4  7.7  92.3
-5  8.2  30.8  69.2
-8  1.2  61.5  38.5
-10  8.2  46.2  53.8

Table 6.3: Error rates for the operator data using δ(x_{n+1}, \overline{x})

These results have several obvious features. Firstly, using \frac{∑ f̂(x_{t})δ(x_{n+1},x_{t})}{∑ f̂(x_{t})} as the outlier detection method is unacceptable due to the failure to detect any signal in all but the easiest cases. Although the type I error rate is very low for most of the SNRs, the failure to detect signals makes this method of very little value.

SNR Type I (%) Type II (%) Sound identified (%)
+5  49.4  0.0  100.0
+1  55.3  0.0  100.0
0  63.5  0.0  100.0
-1  82.4  0.0  100.0
-5  97.6  0.0  100.0
-8  67.1  0.0  100.0
-10  97.6  0.0  100.0

Table 6.4: Error rates for the operator data using δ(x_{n+1}, \overline{x_{w}})

SNR Type I (%) Type II (%) Sound identified (%)
+5  28.2  0.0  100.0
+1  15.3  0.0  100.0
0  5.9  7.7  92.3
-1  14.1  0.0  100.0
-5  14.1  23.1  76.9
-8  5.9  38.5  61.5
-10  10.6  30.8  69.2
Table 6.5: Error rates for the operator data using \overline{δ(x_{n+1}, x_{t})}

Secondly, using δ(x_{n+1}, \overline{x_{w}}) is fraught with problems for entirely the opposite reason. Although the signal detection is the best of any method, it is at the expense of a high type I error rate. Unless it was vital not to miss any part of a signal, it is difficult to see this technique being of general use.

Finally, it is obvious that δ(x_{n+1}, x_{m}) is the best detection method with this test data, combining as it does, minimal type I error rates with a high signal detection result, even at very low signal to noise ratios.

As a comparison, Table 6.8 shows the results when using the Fourier based energy detection method detailed above. It should be noted that although the type I rates are low, the signal detection results are generally much worse than for the kernel methods.

With the actual underwater sounds, it is obviously not possible to calculate type I and type II error rates as the true positions of the sounds are unknown, one way to check results is graphically, by plotting them with representations of the underwater sounds.

SNR Type I (%) Type II (%) Sound identified (%)
+5  25.9  0.0  100.0
+1  7.1  15.4  84.6
0  0.0  23.1  76.9
-1  0.0  100.0  0.0
-5  0.0  100.0  0.0
-8  0.0  100.0  0.0
-10  0.0  100.0  0.0

Table 6.6: Error rates for the operator data using \frac{∑ f̂(x_{t})δ(x_{n+1},x_{t})}{∑ f̂(x_{t})}

SNR Type I (%) Type II (%) Sound identified (%)
+5  17.6  0.0  100.0
+1  11.8  0.0  100.0
0  3.5  7.7  92.3
-1  10.6  0.0  100.0
-5  3.5  0.0  100.0
-8  3.5  0.0  100.0
-10  9.4  15.4  84.6

Table 6.7: Error rates for the operator data using δ(x_{n+1}, x_{m})

Figures 6.6 through 6.9 show the results for the four classes of underwater sound used throughout this work; rain, seals, shrimps and dolphins. These figures show (from the left) levels 11–13 of the wavelet decomposition (these were used to form the summary variables), the results of the six outlier detection methods (kernel value at the new point, squared Mahalanobis distance from the new point to the mean of the points in the density, squared Mahalanobis distance form the new point to the mean of the points in the density weighted by their kernel values, mean squared Mahalanobis distance from the new point to each of the points in the density, squared Mahalanobis distance weighted by the kernel values and the squared Mahalanobis distance from the new point to the point in the density with the maximum kernel value) and the results of the two Fourier based methods (energy detection and spectral detection). All have time on the other axis.

The results of the kernel and Fourier methods are shown as light and dark areas — light areas indicate a signal was detected and dark areas show noise.

Figure 6.6: Results for the analysis of the rain data. 157 Figure 6.7: Results for the analysis of the seal data.

158 Figure 6.8: Results for the analysis of the shrimp data.

159 Figure 6.9: Results for the analysis of the dolphin data.

160

SNR Type I (%) Type II (%) Sound identified (%)
+5  3.5  7.7  92.3
+1  1.2  53.8  46.2
0  0.0  53.8  46.2
-1  4.7  7.7  69.3
-5  1.2  100.0  0.0
-8  1.2  92.3  7.7
-10  3.5  100.0  0.0

Table 6.8: Error rates for the operator data using Fourier based energy detection

As these figures are perhaps not the easiest of results to digest, Tables 6.9 through 6.13 offer a comparison of all the methods described. Values are shown for the percentage of sound identified as signal by each method and some comparisons of the differing detection rates between Fourier and kernel based methods. Table 6.9 shows the percentage of sound that was classified as signal by each of the methods. The most useful results are for the rain data as this is known to contain no signal of interest.

Tables 6.10 through 6.13 all show similar data. The identical classification rate is the percentage of sound that was classified by a kernel based method and a Fourier based method to be the same (either both noise or both signal). Figures are given as a comparison with the energy detection method followed by a comparison with the spectral detection method. Fourier noise and kernel signal gives the percentage of sound that was classified as signal by a kernel based method but as noise by a Fourier based method and, finally, Fourier signal and kernel noise shows the percentage of signal that was classified as noise by a kernel based method but as signal by a Fourier based approach.

Generally, we find that the level of agreement between the best kernel based methods and the Fourier based methods is reasonably high, at around the 70–80% level. Of the remaining sound, the kernel based methods tend to discover more signal than the Fourier methods. There is obviously the possibility that these extra signals are false positive results but on balance we are inclined to the view that the Fourier methods simply have a greater

Method  Rain Seals Shrimps Dolphins
f̂(x_{n+1})  15.9 54.1  37.9  59.6
δ(x_{n+1}, \overline{x})  9.5 54.9  54.3  34.7
δ(x_{n+1}, \overline{x_{w}})  92.8 70.5  74.6  22.2
\overline{δ(x_{n+1}, x_{t})}  24.6 70.5  56.6  74.8
\frac{∑ f̂(x_{t})δ(x_{n+1},x_{t})}{∑ f̂(x_{t})}  0.0 71.6  63.4  72.0
δ(x_{n+1}, x_{m})  19.0 54.9  53.0  68.8
Energy Detection 1.3 62.8  35.4  50.1
Spectral Detection 0.7 54.6  44.9  50.1

Table 6.9: Signal detection rates (all expressed as percentages).

Identical Fourier noise and Fourier signal and
Method  classification  kernel signal  kernel noise
f̂(x_{n+1})  83.8 / 83.8  15.4 / 15.7  0.8 / 0.5
δ(x_{n+1}, \overline{x})  90.2 / 90.2  9.0 / 9.3  0.8 / 0.5
δ(x_{n+1}, \overline{x_{w}})  8.2 / 7.9  91.6 / 92.1  0.2 / 0.0
\overline{δ(x_{n+1}, x_{t})}  75.7 / 75.7  23.8 / 24.1  0.8 / 0.5
\frac{∑ f̂(x_{t})δ(x_{n+1},x_{t})}{∑ f̂(x_{t})} 98.7 / 99.3  0.0 / 0.0  1.3 / 0.7
δ(x_{n+1}, x_{m})  81.0 / 81.0  18.4 / 18.7  0.7 / 0.3

Table 6.10: Comparisons of classification for rain data (expressed as percentages of the total sound). Values are given as energy detection / spectral detection.

Identical Fourier noise and Fourier signal and
Method  classification  kernel signal  kernel noise
f̂(x_{n+1})  87.5 / 84.8  3.4 / 7.4  12.1 / 7.9
δ(x_{n+1}, \overline{x})  84.9 / 84.9  3.6 / 7.7  11.5 / 7.3
δ(x_{n+1}, \overline{x_{w}})  81.1/ 75.6  13.3 / 20.2  5.6 / 4.3
\overline{δ(x_{n+1}, x_{t})}  81.8 / 75.9  13.0 / 0.2  5.2 / 4.1
\frac{∑ f̂(x_{t})δ(x_{n+1},x_{t})}{∑ f̂(x_{t})} 80.3 / 75.1  14.3 / 21.0  5.4 / 3.9
δ(x_{n+1}, x_{m})  84.8 / 83.1  6.2 / 11.1  9.0 / 5.7

Table 6.11: Comparisons of classification for seal data (expressed as percentages of the total sound). Values are given as energy detection / spectral detection.

Identical Fourier noise and Fourier signal and
Method  classification  kernel signal  kernel noise
f̂(x_{n+1})  77.5 / 78.9  12.5 / 7.0  10.0 / 14.1
δ(x_{n+1}, \overline{x})  71.3 / 76.9  23.8 / 16.2  4.9 / 6.9
δ(x_{n+1}, \overline{x_{w}})  44.4 / 50.0  47.4 / 39.8  8.2 / 10.2
\overline{δ(x_{n+1}, x_{t})}  70.0 / 75.6  25.6 / 18.0  4.4 / 6.4
\frac{∑ f̂(x_{t})δ(x_{n+1},x_{t})}{∑ f̂(x_{t})} 64.4 / 70.7  31.8 / 23.9  3.8 / 5.4
δ(x_{n+1}, x_{m})  71.0 / 77.2  23.3 / 15.4  5.7 / 7.4

Table 6.12: Comparisons of classification for shrimp data (expressed as percentages of the total sound). Values are given as energy detection / spectral detection.

Identical Fourier noise and Fourier signal and
Method  classification  kernel signal  kernel noise
f̂(x_{n+1})  78.6 / 79.4  15.4 / 15.0  5.9 / 5.6
δ(x_{n+1}, \overline{x})  57.7 / 57.4  13.3 / 13.4  29.0 / 29.2
δ(x_{n+1}, \overline{x_{w}})  55.8 / 55.9  7.6 / 7.8  36.6 / 36.3
\overline{δ(x_{n+1}, x_{t})}  71.7 / 71.7  26.6 / 26.7  1.7 / 1.7
\frac{∑ f̂(x_{t})δ(x_{n+1},x_{t})}{∑ f̂(x_{t})} 57.1 / 56.8  32.5 / 32.6  10.4 / 10.6
δ(x_{n+1}, x_{m})  74.4 / 74.7  22.2 / 22.1  3.4 / 3.2

Table 6.13: Comparisons of classification for dolphin data (expressed as percentages of the total sound). Values are given as energy detection / spectral detection.

tendency to miss real signal. Our reasons for this view are two-fold. Firstly, the results from the operator data tend to encourage this reading of the situation. In all situations where the signal was not very clear, the Fourier method tended to miss real signal while the best kernel methods, such as the “maximal value” method, determined the true signal with a far higher degree of accuracy. Secondly, the graphical representations of the results does not support the idea that the extra signals are all, in truth, noise. Indeed, for some sounds, such as the shrimp clicks, the extra signals align with visible clicks in the wavelet decomposition with some degree of accuracy.

This general conclusion, that the Fourier based approaches have a higher type II error rate than the best kernel based approaches, is consistent with the results from the operator data set and also with the motivation for moving away from Fourier techniques throughout this work. As noted above, the costs associated with type II errors are much higher than those for type I errors so this reduction in the type II error rate, probably at the expense of some increase in the type I error rate, is generally welcome.

More detailed comparisons of these tables with the results shown, leads to the conclusion that δ(x_{n+1}, \overline{x_{w}}) is the worst method to use as it consistently fails to detect signal (a result borne out by the operator data). Similarly, \frac{∑ f̂(x_{t})δ(x_{n+1},x_{t})}{∑ f̂(x_{t})} also performs relatively poorly, especially with the shrimp and dolphin data. Again, this is consistent with the results from the test data above.

Of the remaining outlier detection methods, the squared Mahalanobis distance from the point with the largest kernel value, δ(x_{n+1}, x_{m}), is the best performer. Of all the methods, it manages to combine the detection of the vast majority of signal flagged by the Fourier methods with extra signal detection in plausible positions. The simple measure of the kernel value at the new point, f̂(x_{n+1}), also performs reasonably although on the shrimp signals especially it is outclassed by the “maximal value” measure.

Of course, it should be remembered that the “maximal value” measure also performed superbly with the operator data as described above. It is therefore less surprising that this performance has been duplicated with real world data sets.

In passing, it is pleasing to discover that the shrimp data set is one where the extra signal detection of the “maximal value” is most noticeable. As one of the reasons for using the wavelet rather than the Fourier transform was its theoretical advantages in capturing the type of short time transients such as shrimp clicks, it is pleasant to see this effect in action.

Finally, mention was made above of the advisability of applying thresholding after the signal identification has taken place. Figures 6.10 through 6.12 demonstrate the effect of this procedure. Three of the underwater sounds analysed above are shown as time series of both the original signal and the final thresholded version with areas of noise removed. The “maximal value” results have been used for this task.

Figure 6.10: The rain data before and after signal identification. Figure 6.11: The seal data before and after signal identification.

Figure 6.12: The shrimp data before and after signal identification. 166

Conclusions and Further Work

7.1 Conclusions

The results presented in this thesis demonstrate that wavelet and kernel based methods can be of use in the analysis of underwater sound data. Chapter 2 argued the case for using the wavelet transform as an alternative to the Fourier transform for moving the sound into the frequency domain for analysis. The wavelet transform is to be preferred over the Fourier transform due to its time localisation properties which allow it to represent signals, especially very short duration transients, in a superior way.

Chapter 4 then showed the effects of different thresholding methods on both some standard test signals and a variety of underwater signals. The standard methods were found to be lacking when used on the underwater signals and this led to a modification being proposed. This modification applied a runs test to the wavelet coefficients to adjust the amount of thresholding that the standard techniques employ.

Results presented show that this modification, especially when used in conjunction with VisuShrink thresholding, produces pleasing audible reconstructions of the underwater sounds whilst simultaneously achieving compression of up to 35% for the data. The results of thresholding the test functions were also satisfactory.

Finally, Chapter 6 presented a method of identifying signals against an evolving background continuum, such as the underwater environment, based on the multivariate density estimation of summaries of the wavelet decomposition. Some possible methods of detecting outliers from this density estimate (which represent signals) were also detailed.

This method was compared with Fourier based techniques, and certain advantages were highlighted. A comparison of the various outlier detection methods was also presented for a known test signal. Finally, identification results were presented for the entire range of underwater signals used in this work and a suggestion as to which outlier detection methods are more suitable for use was made.

The combining of thresholding with the signal identification method was also discussed, with various comparative results presented.

The aim of this work was to consider whether wavelet based techniques are useful in the analysis of underwater sounds. As the results show, there is considerable scope for the use of wavelet based techniques in the identification and noise reduction of such signals — albeit with some modifications to the standard approaches.

We would suggest that RunsThresh is a useful addition to the general choice of thresholds available, although it should be noted that it is only really applicable when used for the analysis of underwater signals or other data with similar wavelet decomposition characteristics.

The use of a multivariate summary of the wavelet decomposition to identify signals is also considered to be a valid approach. Results presented in this work suggest that some features of the underwater sounds, especially short-time transients, are better detected by this method than a Fourier based approach — confirming the superiority of the wavelet transform for capturing such signals.

The method presented, of estimating the probability density of a representation of the wavelet decomposition using a multivariate kernel and then detecting outliers from this estimate, has shown itself to be a useful tool for detecting underwater signals, exceeding the performance of more traditional Fourier based approaches.

7.2 Future Work

There are several obvious areas where there is potential for more analysis in the future. These areas include the choice of summary measures, the choice of wavelet technique and the outlier measure used. These areas will now be discussed in more detail.

Possibly the most obvious area that encourages future study is that of summary measures for the wavelet decompositions. At present, simple summaries of the energy levels at certain scale levels are used. This may not be the most appropriate measure for non-underwater signals or sounds where the background environment is of higher than usual energy. Other measures of the signal may be more useful in these contexts, such as a direct measure of the bandwidth of the sound, various measures of the spectral spread, time periods associated with large and small changes in frequency and amplitude and energy concentration (Fristrup and Walkins 1993). Investigations into whether any of these measures would perform in a superior way to the simple energy measure, either alone or in parallel with other measures, seems a fruitful field for future work.

As part of this whole area, note should also be made of the potential to implement a cross-validation based technique for choosing subsets of these variables before they are analysed. Chapter 6 briefly described why it is preferable to use the smallest subset of variables possible during kernel analysis, involving issues of sample size and speed of calculation.

The next area of interest is the choice of wavelet decomposition technique. This work looked solely at the standard discrete wavelet transform using Daubechies’ wavelets. Other methods that could be used with little change to the entire process presented here include the stationary wavelet transform or wavelet packet transform. Being “best-basis” selection methods, these may enhance the performance, either using the straight decomposition or through summary measures (Saito 1994). Similarly, wavelets other than the Daubechies type used here, such as Chui wavelets, may improve the results.

The area of outlier detection from a multivariate density estimate is also an area that may be improved by more work. Although some of the measures presented have proved to be valuable in the detection of outliers from an estimate of the probability density of the background noise continuum, there is always the possibility that other methods may improve performance still further. The discovery and analysis of such methods would certainly be challenging but should also be a rewarding addition to this underdeveloped area of statistics.

Finally, there is of course the fascinating field of classification, both in its own right and as a final post-process in the general train of events described in this work. Analysing the different possible classification techniques and comparing how they perform as part of this process is a rewarding area which is doubtless worthy of study over a number of years. A brief overview of the area will be given next to encourage further work.

7.2.1 Classification

Although signal detection and noise reduction are the most important areas discussed in this thesis, it would obviously be highly useful if the incoming signals, once discovered, could also be classified as being, for example, a dolphin whistle or a shrimp click.

There are many possible ways this problem could be approached, but a method that combines the two ideas discussed in preceding chapters would have a certain elegance about it. The problem of classification would then simply be the last step in a continual process. This would result in a “black box” approach that would identify the signals in the incoming sound, classify them as belonging to a certain group and finally store them in a compressed format for later review.

This implies that there are three general approaches that could be considered, the main difference being when the compression (by thresholding) should be performed. The choices are firstly, that the initial sound is thresholded before any work is performed; secondly, that the signal is thresholded once it has been identified but before any classification is attempted; and lastly, that the thresholding is only performed once the signal is classified.

Chapter 6 has discussed the issue of thresholding in relation to signal identification. It was argued that thresholding prior to identification would be unwise due to the imperfect nature of existing methods and the consequent possibility of signal loss. However, in relation to classification this problem does not arise. By the classification stage, the signal position is already known and hence, any thresholding would simply alter the signal that is fed to the classifier. As long as a consistent approach is taken, this should not lead to the same degree of problems as for identification.

However, this still does not mean that thresholding prior to classification is a good idea. There is still the potential that thresholding may force two distinct types of signals to appear the same to the classifier, degrading the performance considerably. Conversely, thresholding has the potential to heighten the differences between similar sounds, improving the classification rate. Basically, this means that there is no obvious theoretical justification for or against thresholding prior to classification — only experimentation would show if either method were preferable.

As the end result of the identification process is a multivariate data set, it is natural to consider the more popular multivariate classification techniques that already exist. Such methods include principal component analysis (PCA), discriminant analysis and cluster analysis (Krzanowski 1988).

Another approach that is becoming increasingly popular is that of neural networks, which have already been used in the analysis of underwater sounds as a combined identifier and classifier (Smith et al. 1993). The advantage of pre-identifying signals before using a neural network for classification is that the network is only trying to perform one function, whereas in the case presented in Smith (1993), the network is attempting to both identify and classify the signal.

The last set of possible approaches that will be discussed here can be broadly termed “data dependent”, in that they use natural features of the signal, such as its dominant frequency and length, to classify according to species and type.

Each of these methods will now be briefly examined, with possible approaches being considered in each case.

Multivariate Techniques

The whole area of multivariate discrimination and classification is a wide one that has many texts devoted to it. It would be inappropriate to present here any but the most cursory glance at different possible techniques that could be used. A more detailed analysis of these, and other, methods is given by Krzanowski (1988).

Whilst common techniques such as discriminant analysis, cluster analysis and principal components will doubtless be familiar to many readers, other methods may not. It is this latter class that will be examined here.

Two of the most powerful classification techniques to gain wide acceptance are nearest neighbour and hidden Markov models methods. These will not be detailed in depth but a brief description of both ideas will be presented.

The idea behind the nearest neighbour technique is remarkably simple. A training set of known signals is needed for comparisons (as with most statistical techniques) but if this is obtained then the algorithm is simple.

A distance is calculated from the new observation to all the training observations (using, for example, the squared Mahalanobis distance) and the new observation is assigned to the class which gives the smallest distance measure — the nearest neighbour.

There are problems where multiple sound classes occur simultaneously (Ghosh et al. 1992) and could also fall foul of classes with the odd extreme member. To try and overcome these problems, the k-nearest neighbour method has been developed. In this case, the distances to the training samples are ranked and the k closest ones are examined. The new observation is assigned to the class with most members among the k closest. Such a technique is especially helpful with small training sets and large amounts of noise.

Hidden Markov Models (HMMs) are probably the most successful algorithm used in speech recognition to date, involving as it does an explicit modelling of time (Bahl et al. 1982).

The HMM assumes that the sound is generated by a stochastic system with a finite number of allowable states with certain state transition probabilities — a Markov chain. These transition probabilities govern the evolution of the sound. The most likely series of state transitions to produce the end sound are determined and then, given the model of the stochastic system, the model with the highest probability of generating the observed data is chosen. The model structure must be defined using a priori constraints arising from some understanding of the background process.

Obviously, given the wide range of underwater signals which can be encountered, such model building is a non-trivial task.

Neural Networks

Artificial Neural Networks (ANNs) are an attempt to simulate the process by which the brain makes decisions. In essence, they consist, like the brain, of a set of simple processing units which are linked together with connections of varying weights. In an ANN, these connections are mostly uni-directional, leading to a path through the system from an initial set of inputs to an output. The path is a map from the inputs to the output depending on the weights between units (Ripley 1994).

It is possible that as well as the input and output layers, there are other units in the path which are not accessible from the outside environment. These are known as hidden layers and the justification for them in classification work is given in Smith (1993).

There are obviously many different types of connections that could be made within and between layers, a study of which is beyond the scope of this work. Interested readers are referred to Bishop (1995) which deals with the area of ANNs and pattern recognition.

For the underwater sound problem, the inputs to a neural network could be in one of three forms. Firstly, the raw or thresholded signal could be presented to the network, similar to Smith (1993), but with the advantage that the network would merely be trying to classify rather than identify and classify as given in that work. This should improve the performance somewhat.

Secondly, the multivariate data resulting from the summarisation could be presented to the network as the set of inputs. This would have the advantage of keeping the number of inputs to a reasonable level (in the range 3–8 for most of the underwater sound data) making the network a reasonable size.

Lastly, a “trace vector” for the signal could be presented. This is mearly the multi-dimensional track of the signal from the moment it leaves the kernel estimate and is flagged as an outlier to the moment it returns. For the three dimensional case presented in chapter 6, this would give a set of three dimensional vectors lasting varying lengths of time directly proportional to the signal length. The network would then classify according to these inputs. Again, the number of inputs is kept to a reasonable size.

Data Dependent Methods

Previous chapters of this work have shown that even signals that are similar in the time domain, such as shrimp snaps and dolphin clicks, have different representations when viewed as wavelet decompositions. Other sounds, such as seal and dolphin barks, are differentiated by their bandwidth and duration.

It seems sensible to use some of this prior knowledge as an aid to classification. Such a classifier would obviously be highly dependent on the data being analysed, with an underwater signal classifier having no applicability as a classifier for stress fractures for example. To offset this disadvantage however, it would be readily implementable and would need very little computer power. It is also an obvious backend approach to the “black box” idea advocated above.

Such an approach could be implemented by having a series of checks, one for each “feature” required and then comparing the results of these checks with a known database of signals. There is also the possibility of extracting such important features for each signal and then using these as inputs to either a neural network or classification tree method and classifying in this way.

7.2.2 Discussion

Obviously, these different methods need to be applied to discover which, if any, are most suited to the task of identifying underwater signals, with results balanced against ease of adaption when new signals are required to be classified.

It would appear that the data dependent approach, either implemented alone or in combination with some form of classification tree, offers the potential for very high classification rates together with speed and the possibility of adding new signals, as long as they have unique features when viewed as wavelet decompositions.

Such an approach should exceed the results of the data driven methods as useful a-priori knowledge is brought to bear on the problem. Of course, the other techniques are more flexible when new signals need to be classified.

The final choice of which approach is the most suitable probably depends on the application. If the need to classify new sounds is low compared with the requirement for correct classification of known sounds then some form of data dependent technique would be more suitable.

It should also be noted that given the nature of these techniques, some form of thresholding prior to classification would probably enhance the results. Levelling the field by noise removal should allow similar sounds to be more recognisable by removing extraneous factors that bias the comparison.

Appendix A

Syntax for RunsThresh

RunsThresh is the new method of thresholding the wavelet decomposition detailed in chapter 4. It has been shown to be especially suitable for use with underwater animal sounds such as dolphin whistles and shrimp clicks.

RunsThresh is an adaption of the S-Plus method threshold.wd(), which is part of the WaveThresh 3 package as developed by Nason (1993). The procedure takes the form of a new flag for threshold() — this flag is called

“do.runs”.

A.1 Parameters for RunsThresh

When the flag “do.runs=T” is passed to threshold(), two optional parameters can be used to control the runs test and the result of a positive test (both of these are described in chapter 4).

runs.size

controls the number of coefficients used in one runs test. It is the number of coefficients in the highest scale level of the wavelet decomposition, levels below this have the value progressively halved, down to a minimum of 32 (as explained in chapter 4).

runs.adjust

controls the amount of the original threshold that is used when a test is positive. eġ˙ The default of 0.1 will use 10% of the true threshold.

As can be seen, this formulation allows the runs test to be used with any other threshold method supported by WaveThresh. In particular, the universal policy can be used directly, or in combination with the manual policy to form a threshold value for the entire decomposition based on the highest levels coefficients.

A.2 Example of Usage

This example will show the process of decomposing, RunsThreshing and reconstructing the dolphin whistle shown in Figure 1.1. The wavelet decomposition is shown in Figure 1.2 whilst Figures 1.3 and 1.4 show the decomposition after RunsThresh has been applied and the reconstruction respectively.

The following code should entered into S-Plus after the installation of both WaveThresh 3 and the adapted threshold.wd() function. The whistle should be stored as the object whistle:

> library(WaveThresh)
> whistle.wd <- wd(whistle, family = "DaubLeAsymm",
+           filter.number = 8)
> whT <- threshold(whistle.wd, type = "hard",
+           policy = "universal", dev = mad.est,
+            return.threshold = T,
+           levels = nlevels(whistle.wd)-1)
> whistle.run <- threshold(whistle.wd, type = "hard",
+                   policy = "manual", levels =
+                   start.level:(nlevels(whistle.wd) - 1),
+                   value = rep(whT, nlevels(whistle.wd)
+                   - start.level), do.runs = T,
+                   runs.adjust = 0.1, runs.size = 128)

The first command simply starts WaveThresh 3 whilst the second command performs the DWT. The family and filter.number allow control over the wavelet used for this decomposition. The next two commands actually perform the thresholding. The first finds the VisuShrink threshold based on the lowest level of coefficients to be used as the basis for RunsThresh. The mad.est() function simply calculates the adjusted median absolute deviation suggested in Donoho and Johnstone (1994).

The final command actually performs the RunsThresh command. The important parameters are the last three which specify a runs test of size 128 (on the whistle, which is 32768 data points, this gives 128 windows at the highest scale level) and an adjustment to 10% of the VisuShrink threshold value for coefficient groups with a positive runs test.

A.3 Credits

This work is heavily based on the WaveThresh 3 code described in Nason and Silverman (1994) which was most beneficial. The runs test procedure was detailed in Siegel (1956).

Appendix B

Syntax for KP

KP — Kernel Program — is used to implement the ideas presented in Chapter 6. It is C++ code to analyse a matrix of summary values (consisting of known noise and new observations) and determine if the new observations are noise or signal.

This program was used to obtain the results given in Chapter 6 for the signal identification using the kernel method. It works by forming a probability density estimate on the initial columns of the input matrix, using a Fukunaga (1972) kernel, and then comparing each new observation (column of the input matrix) with the kernel to discover if it is also noise using the methods described below.

B.1 Starting The Program

The syntax for the program is straight forward, most of the options, which are described below, being set after the start. At the command line type:

   kp datafile resultfile

datafile is the standard matrix file as described below and resultfile is where the results are to be placed.

# Optional comment lines starting with a hash rows = x columns = y # The above can be written as r = x c = y ::: 1 2 3 2 3 4

Figure B.1: The standard matrix file format

B.1.1 The Standard Matrix File Format

The standard matrix file format is defined in Eckel (1993) for use with his matrix manipulation routines. The basic format of the file is shown in figure B.1 where the data (given after the ::: line) are presented row-wise. This format is easy to create — S-Plus for instance can create such files with ease.

B.1.2 The Results File Format

The results file output by the program has a simple matrix format. Each observation has eight characteristics given on one row of the output file. These are the observation number, the five values for the distance measured described in Chapter 6, the value that the distance measure is compared with and the result (given as 1 for signal and 0 for noise).

B.2 Parameters Used By KP

On starting the program, a list of parameters are obtained from the operator. These determine various modes of operation and analysis. A detailed list of these options will now be given, in the order that they are requested.

How much junk at start of file?

This allows number of columns at the beginning of the input matrix to be completely ignored.

How much data for kernel?

This number should be the number of columns of the input matrix which are to be used for the initial noise segment — it is also the number of observations used in the kernel calculations.

How many rows of data should be used?

This determines the number of variables to be used in the calculations — it must not exceed the number of rows of the input matrix.

Type in row numbers separated by spaces.

Here, the actual rows of the input matrix that contain the variables are given. The first row is represented by 1, the second by 2, etc.

How many signals should be stored?

This determines how many signals are held for later display. At present, this option does nothing.

What size buffer?

The size of the buffer to use with 0 representing no buffer. This is fully explained in Chapter 6.

What size bandwidth?

Allows a fixed bandwidth to be used, if 0 is entered then a crossvalidated bandwidth is used, as described in Chapter 6.

Outlier measure.

The method for determining outliers is specified here. Details of the various methods used are given in Chapter 6. The possible responses are k for the kernel value, d for the distance, w for the weighted distance, m for the mode (not yet implemented), n for the mean and t for the weighted mean.

What percentage cutoff?

This determines how extreme the observation must be to be considered an outlier. The description in Chapter 6 for the outlier detection methods uses minimum or maximum values from the existing noise as comparisons but this parameter allows a value to be chosen which is a percentage into the ordered list of values.

Update bandwidth every x observations.

How often should the bandwidth be updated. In theory, a new bandwidth should be calculated every observation, in practice however, computational speed often dictates that some gap should occur before the time intensive bandwidth calculations are performed. This option only applies if cross-validated bandwidths are used.

How many input files?

Multiple input files can be used, to aid the analysis of continuous data streams. The only restriction is that they must all have the format basename.x where x increases from 1.

Multiply all input by?

In some applications, the input variables may all be very large or small. In such circumstances, it may be helpful to pre-scale all observations by some fixed factor. This option allows such a process.

B.3 Problems With Implementation

The actual method used for this analysis is very promising, Chapter 6 shows some comparative results, but the actual implementation has some problems.

The major problem is the speed. The program is much too slow — 12 seconds takes approximately six weeks to analyse in full on a SPARCclassic with 16 MB of RAM. There are two possible reasons for this, the first is that the kernel used is a normal density. As pointed out in Silverman (1986), this is “foolish” as it requires a call to the exponential function for each evaluation. As the bandwidth is calculated by cross-validation under the automatic scheme, many such evaluations are required each iteration.

The second possible reason for the slowness is a general lack of optimisations in the code. With more attention being paid to such coding issues, the program could be made considerably quicker. As a minor point, the actual computer used to run the program on is very slow and a faster computer would also be a great help.

As noted by Silverman (1986), the speed would be greatly enhanced if the kernel were calculated using a fast Fourier transform based method as detailed in Silverman (1982). This is probably the most appropriate method of enhancing the speed of this program.

In actual fact, experimentation has shown that the main actual causes of the slowness of the program are the value of ϵ used in the golden section search for the optimal bandwidth and updating the bandwidth for every new observation. When ϵ is reduced to a value of, for instance, 0.00001 and the bandwidth is only updated every 25 new observations the program completes the analysis in 50 minutes. As the actual bandwidths used with the underwater sounds tend to be approximately 3 and only fluctuate slowly over time, this is a reasonable compromise. It would however be useful if ϵ values could be set interactively, as with the other parameters.

As this single change produces results 1000 times faster than the original method, it was considered unnecessary to make any of the other possible enhancements as they would not produce such dramatic improvements.

B.4 Credits

Silverman (1986) and Wand and Jones (1995) provided the description of the kernel method for probability density estimation. Tukey and Tukey (1981) gave details of “best-practice” covariance matrix implementation, example code was provided by Herraman (1968) and Clarke (1971). The C++ matrix class is from Eckel (1993) and is used with permission.

B.5 The Program

The program was complied and run on a SPARCclassic with 16 Mb of RAM using the GNU C++ complier. It is copyright 1996-7 by Ken Powell but may be used for any non-commercial purpose. Other users please contact me first.

The code is available from the ftp site msor.ex.ac.uk in the directory /pub/kpowell with the filename kernel.tgz. It compiles with gcc version 2.7.2 or later but needs the matrix class code from Eckel (1993) which is copyright of the book’s author.

Bibliography

Abramovich, F. and Y. Benjamini (1995). Thresholding of wavelet coefficients as multiple hypothesis testing procedure. In A. Antoniadas and G. Oppenheim (Eds.), Wavelets and Statistics, Volume 103 of Lecture Notes in Statistics, pp. 5–14. New York: Springer-Verlag.

Abramovich, F. and Y. Benjamini (1996). Adaptive thresholding of wavelet coefficients. Computational Statistics & Data Analysis 22, 351–361.

Abramovich, F., T. Sapatinas, and B. W. Silverman (1998). Wavelet thresholding via a bayesian approach. Journal of the Royal Statistical Society (Series B) 60, 000–000.

Abramovich, F. and B. W. Silverman (1998). Wavelet decomposition approach to statistical inverse problems. Biometrika 85, 000–000.

Altman, N. S. (1990). Kernel smoothing of data with correlated errors. Journal of the American Statistical Association 85 (411), 749–759.

Antoniadis, A. and G. Oppenheim (1995). Wavelets and Statistics, Volume 103 of Lecture Notes in Statistics. New York: Springer-Verlag.

Atkinson, A. C. and H.-M. Mulira (1993). The stalactite plot for the detection of multivariate outliers. Statistics and Computing 3, 27–35.

Atlas, L., R. Cole, Y. Muthusamy, Y. Lippman, J. Conner, D. Park, M. El-Sharkawi, and R. J. M. II (1990). A performance comparison of trained multilayer perceptrons and trained classification trees. In C. Lau (Ed.), Neural Networks; Theoretical Foundations and Analysis. IEEE Press.

Auty, L. A. (1986). Study of the whales and dolphins occurring in the north atlantic: Including distribution, migration and sound production. BSc nautical studies honours project, Polytechnic South West, Plymouth.

Bahl, L. R., F. Jelinek, and R. L. Mercer (1982). Continuous speech recognition: Statistical methods. In L. N. Kanal and P. R. Krishnaiah (Eds.), Handbook of Statistics Vol. 2, pp. 549–573. North Holland Publishing Company.

Bailey, T. C., T. Sapatinas, K. J. Powell, and W. J. Krzanowski (1996). Use of wavelets in the detection of underwater sound signals. In Interface96, Sydney, Australia.

Bailey, T. C., T. Sapatinas, K. J. Powell, and W. J. Krzanowski (1998). Signal detection in underwater sound using wavelets. Journal of the American Statistical Association 93, 000–000. To appear March 1998.

Benedetto, J. J. and M. W. Frazier (1994). Wavelets: Mathematics and Applications. Studies in Advanced Mathematics. Boca Raton: CRC Press.

Bishop, C. M. (1995). Neural Networks for Pattern Recognition. Oxford: Claredon Press.

Boashash, B. (1990). Time-frequency signal analysis. In S. Haykin (Ed.), Advances in Spectrum Analysis and Array Processing, pp. 418–517. Englewood Cliffs, New Jersey: Prentice-Hall.

Boashash, B. and P. O’Shea (1990). A methodology for detection and classification of some underwater acoustic signals using time-frequency analysis techniques. IEEE Transactions on Acoustic, Speech and Signal Processing 38, 1829–1841.

Bock, M. E. (1992). Estimating functions with wavelets. Statistical Computing and Statistical Graphics Newsletter 3 (1), 4–8.

Bock, M. E. and G. Pliego (1992). Estimating functions with wavelets. part II: Using a daubechies wavelet in nonparametric regression. Statistical Computing and Statistical Graphics Newsletter 3 (2), 27–34.

Boutell, T. (1996, October). PNG (portable network graphics) specification - W3C recommendation. Technical report, PNG Development Group, WWW Consortium.

Bowman, A. W. (1984). An alternative method of cross-validation for the smoothing of density estimates. Biometrika 71, 353–360.

Bruce, A. and H.-Y. Gao (1994). S+WAVELETS User’s Manual, Version 1.0. Seattle: StatSci, a division of MathSoft, Inc.

Chen, C. H. (1985). Recognition of underwater transient patterns. Pattern Recognition 18, 485–490.

Chipman, H. A., E. D. Kolaczyk, and R. E. McCulloch (1996). Adaptive bayesian wavelet shrinkage. Technical Report 421, Dept. of Statistics, Chicago University, USA.

Chui, C. K. (1992). An Introduction to Wavelets. New York: Academic Press.

Clarke, M. R. B. (1971). AS 41: Updating the sample mean and dispersion matrix. Appl. Statist. 20 (2), 206–209.

Clyde, M., G. Parmigiani, and B. Vidakovic (1995). Multiple shrinkage and subset selection in wavelets. Discussion Paper 95-37, ISDS, Duke University, USA.

Cohen, L. (1989). Time-frequency distributions — a review. Proceedings of the IEEE 77, 941–981.

Coifman, R. R. and N. Saito (1994). Construction of local orthonormal bases for classification and regression. Comptes Rendus des séances de l’Acadèmie Sciences, Série I 319, 191–196.

Cooley, J. W. and J. W. Tukey (1965). An algorithm for the machine computation of complex fourier series. Meth Comput 19, 297–301.

Creasey, D. J. (1985). Advanced Signal Processing. IEE Telecommunications Series. Exeter: The Institution of Electrical Engineers.

Creasey, D. J., B. V. Smith, and B. K. Gazey (1989a). Underwater acoustics and sonar systems - acoustics, transducers and propagation. Lecture Notes Vol. 1, Birmingham University.

Creasey, D. J., B. V. Smith, and B. K. Gazey (1989b). Underwater acoustics and sonar systems - signal processing and systems. Lecture Notes Vol. 2, Birmingham University.

Daubechies, I. (1990). The wavelet transform, time-frequency localization and signal analysis. IEEE Transactions on Information Theory 36 (5), 961–1005.

Daubechies, I. (1992). Ten Lectures on Wavelets. Philadelphia: Society for Industrial and Applied Mathematics.

Daubechies, I. (Ed.) (1993). Different Perspectives on Wavelets, San Antonio, Texas. American Mathematical Society, Rhode Island.

Donoho, D. L. (1993a). Nonlinear wavelet methods for recovery of signals, densities, and spectra from indirect and noisy data. In I. Daubechies (Ed.), Proceedings of Symposia in Applied Mathematics, Volume 47.

Donoho, D. L. (1993b, January). Wavelet shrinkage and W.V.D.: A 10minute tour. Technical Report 416, Dept. of Statistics, Stanford University, Stanford, California.

Donoho, D. L. and I. M. Johnstone (1994). Ideal spatial adaption by wavelet shrinkage. Biometrika 81 (3), 425–455.

Donoho, D. L. and I. M. Johnstone (1995). Adapting to unknown smoothness via wavelet shrinkage. Journal of the American Statistical Association 90 (432), 1200–1224.

Donoho, D. L. and I. M. Johnstone (1997). Minimax estimation via wavelet shrinkage. Annals of Statistics 25, 000–000. (In press).

Donoho, D. L., I. M. Johnstone, G. Kerkyacharian, and D. Picard (1995). Wavelet shrinkage: Asymptopia? (with discussion). Journal of the Royal Statistical Society (Series B) 57 (2), 301–369.

Drumheller, D. M. (1991, July). Theory and application of the wavelet transform to signal processing. Technical Report 9316, Naval Research Laboratory, Washington, DC.

Duin, R. P. W. (1976). On the choice of smoothing parameters for Parzen estimators of probability density functions. IEEE Transactions on Computers C-25, 1175–1179.

Eckel, B. (1993). C++ Inside and Out. Berkley, California: Osborne McGraw-Hill.

Friedman, J. H., F. Baskett, and L. J. Shustek (1975). An algorithm for finding nearest neighbours. IEEE Transactions on Computers 24, 1000–1006.

Fristrup, K. M. and W. A. Walkins (1993, October). Marine animal sound classification. Technical report, Woods Hole Oceanographic Institution.

Fuchs, C. and Y. Benjamini (1994). Multivariate profile charts for statistical process control. Technometrics 36, 182–195.

Fukunaga, K. (1972). Introduction to Statistical Pattern Recognition. New York: Academic Press.

Fulton, J., G. Renda, and D. Lemke (1994). Overview of NCDaudio. Network Computing Devices, Inc.

Ghosh, J., L. M. Deuser, and S. D. Beck (1992). A neural network based hybrid system for detection, characterization and classification of shortduration oceanic signals. IEEE Journal of Oceanic Engineering 17, 351–363.

Gorman, R. P. and T. J. Sejnowski (1988). Analysis of hidden units in a layered network trained to classify sonar targets. Neural Networks 1, 75–89.

Grant, P. M., C. F. N. Cowan, B. Mulgrew, and J. H. Dripps (1989). Analogue and Digital Signal Processing and Coding. Sweden: Chartwell-Bratt.

Habbema, J. D. F., J. Hermans, and K. van der Brock (1974). A stepwise discrimination program using density estimation. In G. Bruchman (Ed.), Compstat 1974, pp. 100–110. Vienna: Physica Verlag.

Hall, P. and G. P. Nason (1997). On choosing a non-integer resolution level when using wavelet methods. Statistics & Probability Letters 34, 5–11.

Hart, J. D. (1991). Kernel regression estimation with time series errors. Journal of the Royal Statistical Society (Series B) 53 (1), 173–187.

Hermans, J., J. D. F. Habbema, T. K. D. Kasanmoentalib, and J. W. Raatgever (1982, March). Manual for the Alloc80 Discriminant Analysis Program. Leiden, Netherlands: Dept. of Medical Statistics, University of Leiden.

Herraman, C. (1968). AS 12: Sums of squares and products matrix. Appl. Statist. 17 (3), 289–292.

Hopper, T. and F. Preston (1992). Compression of grey-scale fingerprint images. In Data Compression Conference, New York. IEEE Computer Society Press.

Johnstone, I. M. and B. W. Silverman (1997). Wavelet threshold estimators for data with correlated noise. Journal of the Royal Statistical Society (Series B) 59 (2), 319–351.

Jones, D. L. and T. W. Parks (1992). A resolution comparison of several time-frequency representations. IEEE Transactions on Signal Processing 40, 413–420.

Jones, M. C. and H. W. Lotwick (1984). A remark on algorithm as 176. kernel density estimation using the fast fourier transform. Applied Statistics 33 (1), 120–122. AS R50.

Kay, J. and B. Kummerfeld (1988). C Programming in a UNIX Environment. International Computer Science Series. Sydney: Addison-Wesley Publishing Company.

Khattree, R. and D. N. Naik (1987). Detection of outliers in bivariate time series data. Communications in Statistics — Theory and Methods 16, 3701–3714.

Kidd, R. C. (1995). Comparison of wavelet scalar quantization and jpeg for fingerprint image compression. Journal of Electronic Imaging 4 (1), 31–39.

Krzanowski, W. J. (1988). Principles of Multivariate Analysis - A User’s Perspective. Number 3 in Oxford Statistical Science Series. Oxford: Oxford University Press.

Krzanowski, W. J. and T. C. Bailey (1996, May). Multivariate statistical pre-processing of underwater sound data. Agreement No. 2041/045 for DRA, Winfrith Technology Centre, Dorset.

Lammi, L. and T. Sarjakoski (1995). Image compression by the JPEG algorithm. Photogrammetric Engineering and Remote Sensing 61 (10), 1261–1266.

Learned, R. E. and A. S. Wilsky (1995). A wavelet packet approach to transient signal classification. Applied and Computational Harmonic Analysis 2, 265–278.

Lippman, S. B. (1991). C++ Primer (second ed.). Reading, Massachusetts: Addison-Wesley Publishing Company.

Ludeman, L. C. (1986). Fundamentals of Digital Signal Processing. Ellis Horwood series in Computers and their Applications. Cambridge, USA: Harper and Row.

Mallat, S. (1989). A theory for multiresolution signal decomposition: the wavelet representation. IEEE Transactions on Pattern Analysis and Machine Intelligence 11 (7), 674–693.

Marshall, A. G. and F. R. Verdun (1990). Fourier Transforms in NMR, Optical, and Mass Spectrometry: A User’s Handbook. Amsterdam: Elsevier.

McCausland, I. (1969). Introduction to Optimal Control. New York: John Wiley and Sons, Inc.

Meyer, Y. (1993a). Book review: An introduction to wavelets and ten lectures on wavelets. Bulletin of the American Mathematical Society 28 (2), 350–360.

Meyer, Y. (1993b). Wavelets: Algorithms and Applications. Philadelphia: Society for Industrial and Applied Mathematics. Translated by Robert D Ryan.

Nason, G. P. (1993). The WaveThresh package; wavelet transform and thresholding software for S. Available from the StatLib archive.

Nason, G. P. (1996). Wavelet shrinkage using cross-validation. Journal of the Royal Statistical Society (Series B) 58, 463–479.

Nason, G. P. (1997). Wavelets and statistics: Wavethresh 3. Unpublished book.

Nason, G. P. and B. W. Silverman (1994). The discrete wavelet transform in S. Journal of Computational and Graphical Statistics 3, 163–191.

Nason, G. P. and B. W. Silverman (1995). The stationary wavelet transform and some statistical applications. In A. Antoniadis and G. Oppenheim (Eds.), Wavelets and Statistics, Volume 103 of Lecture Notes in Statistics, pp. 281–300. New York: Springer-Verlag.

Nussbaumer, H. J. (1981). Fast Fourier Transform and Convolution Algorithms. Springer Series in Information Sciences. Berlin: Springer-Verlag.

Ogden, T. and E. Parzen (1996a). Change point approach to data analytic wavelet thresholding. Statistics and Computing 6, 93–99.

Ogden, T. and E. Parzen (1996b). Data dependent wavelet thresholding in nonparametric regression with change point applications. Computational Statistics & Data Analysis 22, 53–70.

Pal, N. R. and S. K. Pal (1993). A review on image segmentation techniques. Pattern Recognition 26, 1277–1294.

Park, B. U. and J. S. Marron (1990). Comparison of data-driven bandwidth selectors. Journal of the American Statistical Association 85, 66–72.

Park, B. U. and B. Turlach (1992). Performance of several data driven bandwidth selectors (with discussion). Computational Statistics 7, 251– 285.

Penny, K. I. (1996). Appropriate critical values when testing for a single multivariate outlier by using the mahalanobis distance. Applied Statistics 45, 73–81.

Pesquest, J. C., H. Krim, and H. Corfanban (1994). Time invariant orthonormal wavelet representations. Submitted for publication.

Pettit, F. (1985). Fourier Transforms in Action. Bromley: Chartwell-Bratt.

Powell, K. J., T. Sapatinas, T. C. Bailey, and W. J. Krzanowski (1995). Application of wavelets to the pre-processing of underwater sounds. Statistics and Computing 5, 265–273.

Prasad, S., M. Chakraborty, and H. Parhasarathy (1995). The role of statistics in signal processing — a brief review and some emerging trends. Indian Journal of Pure and Applied Mathematics 26, 547–578.

Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery (1992). Numerical Recipes in C (second ed.). Cambridge: Cambridge University Press.

Remme, J., J. D. F. Habbema, and J. Hermans (1980). A simulative comparison of linear, quadratic and kernel discrimination. Journal of Statistical Computation and Simulation 11, 87–106.

Ripley, B. D. (1994). Neural networks and related methods for classification. Journal of the Royal Statistical Society (Series B) 56 (1), 409–456.

Saito, N. (1994, December). Local Feature Extraction and its Applications Using a Library of Bases. Ph. D. thesis, Yale University, Yale, USA.

Scott, D. W. (1992). Multivariate Density Estimation: Theory, Practice and Visualization. Wiley Series in Probability and Mathematical Statistics. New York: John Wiley and Sons, Inc.

Siegel, S. (1956). Nonparametric Statistics for the Behavioural Sciences. Tokyo: McGraw-Hill Kogakusha.

Silverman, B. W. (1982). Kernel density estimation using the fast fourier transform. Applied Statistics 31 (1), 93–99. AS 176.

Silverman, B. W. (1986). Density Estimation for Statistics and Data Analysis. Number 26 in Monographs on Statistics and Applied Probability. London: Chapman and Hall.

Smith, D. J. (1993, April). The Application of Artificial Neural Networks to the Classification of Underwater Cetacean Sounds. Ph. D. thesis, Exeter University, Exeter.

Smith, D. J., T. C. Bailey, and A. G. Munford (1993). Robust classification of high dimensional data using artificial neural networks. Statistics and Computing 3, 71–81.

Smith, M. J. T. and R. M. Mersereau (1992). Introduction to Digital Signal Processing. New York: John Wiley and Sons, Inc.

Statistical Sciences, Inc (1993a). S-PLUS Guide to Statistical and Mathematical Analysis, Version 3.2. Seattle: StatSci, a division of MathSoft, Inc.

Statistical Sciences, Inc (1993b). S-PLUS Programmer’s Manual, Version 3.2. Seattle: StatSci, a division of MathSoft, Inc.

Statistical Sciences, Inc (1993c). S-PLUS Reference Manual, Version 3.2. Seattle: StatSci, a division of MathSoft, Inc.

Statistical Sciences, Inc (1993d). S-PLUS User’s Manual, Version 3.2. Seattle: StatSci, a division of MathSoft, Inc.

Stein, C. (1981). Estimation of the mean of a multivariate normal distribution. Annals of Statistics 9, 1135–1151.

Strang, G. (1989). Wavelets and dilation equations: A brief introduction. Society for Industrial and Applied Mathematics Review 31 (4), 614–627.

Strang, G. (1993). Wavelets transforms versus fourier transforms. Bulletin of the American Mathematical Society 28 (2), 288–304.

Telfer, B. A., H. H. Szu, G. J. Garcia, H. Ko, A. Dubey, and N. Witherspoon (1994). Adaptive wavelet classification of acoustic and backscatter and imagery. Optical Engineering 33, 2192–2203.

Tizzard, K. (1986). C for Professional Programmers. Ellis Horwood series in Computers and their Applications. Chichester: Ellis Horwood.

Tukey, P. A. and J. W. Tukey (1981). Graphical display of data sets in 3 or more dimensions. In V. Barnett (Ed.), Interpreting Multivariate Data, pp. 189–275. Chichester: John Wiley and Sons.

Urick, R. J. (1975). Principles of Underwater Sound. McGraw Hill.

Vidakovic, B. (1994). Non-linear wavelet shrinkage with bayes rules and bayes factors. Discussion Paper 94-24, ISDS, Duke University, USA.

von Sachs, R., G. P. Nason, and G. Kroisandt (1996). Spectral representation and estimation for locally stationary wavelet processes. In Proceedings of the workshop: Spline functions and wavelets, Montreal.

Walker, J. S. (1988). Fourier Analysis. Oxford: Oxford University Press.

Walker, J. S. (1991). Fast Fourier Transforms. Studies in Advanced Mathematics. Baco Raton: CRC Press.

Wand, M. P. and M. C. Jones (1994). Multivariate plug-in bandwidth selection. Computational Statistics 9, 97–116.

Wand, M. P. and M. C. Jones (1995). Kernel Smoothing. Number 60 in Monographs on Statistics and Applied Probability. London: Chapman and Hall.

Wang, Y. (1996). Function estimation via wavelet shrinkage for long memory data for data. Annals of Statistics 24, 466–482.

Weyrich, N. and G. T. Warhola (1994). De-noising using wavelets and cross-validation. Technical Report AFIT/ENC/TR/94-01, Dept. of Mathematics and Statistics, Air Force Institute of Technology, AFIT/ENC, 2950 P ST, Wright-Patterson Air Force Base, Ohio, 45433-7765.

Whalen, A. D. (1971). Detection of Signals in Noise. Electrical Series. London: Academic Press.

Wiener, R. S. and L. J. Pinson (1988). An Introduction to Object-Orientated Programming and C++. Reading, Massachusetts: Addison-Wesley Publishing Company.

Wierda, S. J. (1994). Multivariate statistical process control — recent results and directions for future research. Statistica Neerlandica 48, 147–168.