Classification of sperm whale clicks (Physeter macrocephalus) with Gaussian-kernel based networks

Algorithms 2009,2, 1232 - 1247; doi:10.3390/a2031232 OPEN ACCESS algorithms ISSN 1999-4893 www.mdpi.com/journal/algorithms Article Classification of Sperm Whale Clicks (Physeter Macrocephalus) with Gaussian-Kernel-Based Networks Mike van der Schaar 1,?, Eric Delory 2and Michel Andr´ e1,? 1Laboratori d’Aplicacions Bioac´ ustiques, Universitat Polit` ecnica de Catalunya, Rambla Exposici´ o s/n, 08800 Vilanova i la Geltr´ u, Spain 2ETIS, UMR 8051 (CNRS, ENSEA, UCP), Avenue du Ponceau 6, BP 44, F-95014 Cergy-Pontoise Cedex, France; E-Mail: [email protected] ?Author to whom correspondence should be addressed; E-Mails: mike.v[email protected] (M.V.); [email protected] (M.A.); Tel.: (34) 938967227; Fax: (34) 938967201. Received: 9 July 2009; in revised form: 31 August 2009 / Accepted: 15 September 2009 / Published: 22 September 2009 Abstract: With the aim of classifying sperm whales, this report compares two methods that can use Gaussian functions, a radial basis function network, and support vector machines which were trained with two different approaches known as C-SVM and ν-SVM. The methods were tested on data recordings from seven different male sperm whales, six containing single click trains and the seventh containing a complete dive. Both types of classifiers could distinguish between the clicks of the seven different whales, but the SVM seemed to have better generalisation towards unknown data, at the cost of needing more information and slower performance. Keywords: classification; sperm whale; radial basis function; support vector machine 1. Introduction Sperm whales (Physeter macrocephalus), when living in a social community, often forage in small groups. During their feeding dive, which may be to depths up to 2 kilometres [1], they start producing sonar signals fairly soon after the start of a dive and generally continue until the ascent back to the surface. Usually one click per second is produced on average, but at times this frequency is increased (presumably when they have found prey) and up to fifty signals per second may be produced [2]. This Algorithms 2009,21233 sequence of very rapid clicks is called a creak, and the period from the start of a click sequence to a creak, or a prolonged moment of silence, is called a click train. The recording of these diving groups results in a mixture of signals, and the manual assignment of a click from a click train to the animal that produced it is often a difficult and arduous task. To this end, we want to build an automated method that can distinguish between clicks from different animals using primarily characteristic information in the clicks themselves. This could then be combined with, for example, time delays of arrival at the hydrophones to reliably reconstruct the original click trains for the individual whales. When recorded on axis, a click consist of a series of up to five pulses separated by a few milliseconds where the second pulse dominates strongly over the others. A relative difference of 40 dB in signal level can be encountered [3]. The clicks are broad band with energy in frequencies exceeding 30 kHz and a peak frequency around 12 kHz. However, when recording off axis, as is almost always the case in practice, the signal becomes distorted with strong attenuation of this second pulse. The click’s characteristics may change depending on the animal’s orientation and distance to the hydrophone, while the animal’s depth [4] and activities may play a role as well. A typical example of a low-pass filtered clicks is shown in Figure 1. Each image in this figure contains 10 consecutive clicks superimposed on each other. Synchronisation was done with a low frequency matched filter. It can be seen how the signal changes through the click train. Focusing on for example the area around sample 210, a new peak appears probably due to different time delays between the click’s pulses that can be caused by changes in the animal’s orientation with respect to the hydrophone. Research has shown [5] that dominant frequencies of off axis clicks for both male and female sperm whales can be found below 2,000 Hz. Since low frequencies are less influenced by orientation or distance, they are more suitable to be searched for constant characteristics. Figure 1. Example clicks from one animal at different moments in the click train. Each image contains 10 consecutive superimposed clicks that were filtered below 3,000 Hz and normalised in energy. A pulse seems to be moving from left to right, especially visible around sample 210. 100 200 300 400 500 −0.2 0 0.2 100 200 300 400 500 −0.2 0 0.2 Samples 100 200 300 400 500 −0.2 0 0.2 An earlier attempt has been made to identify a whale based on modelling these dominant frequencies directly using a Gabor function [6]. The use of Gabor functions is interesting as it suggests that nature uses a signal optimised in a time-frequency sense, and they have been used successfully to describe dolphin sonar [7,8]. However, it was found that the dominant frequencies were not stable enough to be used as an identifier for the entire duration of a dive. Similarly, in the same study, a simple linear classification method was not able to distinguish clearly between individuals in a small group of whales, and therefore this report looks at non-linear classifiers in the form of a neural network. Neural networks Algorithms 2009,21234 have been used in the past for marine mammal classification with some success [9,10]. The distribution of the characteristics from the sperm whale clicks suggested the use of a Gaussian model, and therefore a radial basis function network (RBF) architecture was used in [11] to separate sperm whales. Here, we compare the performance of radial basis functions with support vector machines (SVM) [12,13]. SVM use a similar network architecture, but follow a different underlying approach and are trained differently. Training of the SVM was done in the standard approach, known as C-SVM, and a second approach known as ν-SVM [14,15]. An advantage of the latter method is that there is a more direct control on the number of support vectors used by the machine. 2. Data Acquisition, Preparation and Feature Selection 2.1. Data acquisition The sperm whale data were collected from an inflatable boat during four field seasons spanning four to ten weeks each (from 1997 to 1999) at Kaikoura, New Zealand [16]. Recordings were made of solitary diving male sperm whales using an omni-directional hydrophone (Sonatech 8,178; frequency response 100 Hz to 30 kHz ±5 dB) lowered to a depth of 20 m. This hydrophone was first connected to a fixed gain amplifier (flat response from 0 to 45 kHz) and then to one channel of a Sony TCD-D10PROII Digital Audio Tape recorder (frequency response 20 Hz to 22 kHz ±1 dB with an anti-alias filter at 22 kHz). The recordings were digitized at 48 kHz and 16 bits. The use of data from solitary diving whales guarantees that no data from different animals were mixed, thus helping to obtain optimal results. The typical duration of the recorded click trains was around 2.5 mins. The 30 mins complete dive was a sequence of such click trains. The dive was divided in 10 data segments to ease data handling and manual analysis, but these did not exactly cover 10 click trains as the trains themselves are not always very well defined. A pause in the click production is sometimes too short to consider the continuation to be a new click train, but rather it can indicate that a few signals were not detected. Generally, a click train ends with a creak where it is assumed the animal is capturing a prey, but this is not always heard on a recording as the signal can be too weak. 2.2. Data preparation and feature selection In preparation for the classification algorithm, the clicks were manually detected, filtered for echoes, and checked for acceptable noise levels. These clicks were then denoised using a standard soft-thresholding algorithm, available in Wavelab [17], and synchronised using a matched filter on the low dominant frequency with a typical example click. Initially, the data were band-pass filtered between 100 and 20,000 Hz. Data from seven different animals were available for this study, comprised of six single click trains and a complete dive. This dive was considered to be especially interesting as it allowed to see the performance of the algorithm, and validity of features, for the duration of an entire dive. Therefore, the dive was split up in two unequal parts. One click train early in the dive was separated and joined together with the other six available click trains used for training. The remainder of the dive was put in its own set and was only used to test the classifier, it was never used for training or parameter selection. This approach simulated the situation where a classifier would have to be trained with data at the start of a Algorithms 2009,21235 recording and allowed to assess its capacity to generalise to patterns much later in the dive sequence, that may have undergone changes as in Figure 1. The seven click trains were used to train the classifiers. As the objective of the classification is that a classifier can be trained with only the start of a recording and then autonomously classify the remainder, only the first 50 clicks were used from the start of each click train, which corresponds to roughly 50 seconds. The other available patterns in the click trains were used for validation, but were never considered for training. The features were selected using a local discriminant basis [18,19] for the seven classes. The exact same procedure was followed as in [6]. First, each click in the training set was expressed in a wavelet packet table. Where the usual wavelet filter retains the high-pass wavelet coefficients and continues filtering the low-pass scale coefficients, the wavelet packet table filters both outputs again, creating a redundant library of bases that can be selected for reconstruction of the signal (a detailed discussion about the relationship between wavelets and filter banks can be found in [20]). With each pass through the wavelet filter the frequency band of the input signal is split in two, producing low and high frequency outputs. These outputs are stored and passed through the filter again. In this paper these recursive steps will be called the splitting level, e.g., at level 3 the signal has been passed through the filter twice. At a splitting level lthere will be 2l−1band limited signals, each with a bandwidth of Fs/2lwith Fsthe sampling frequency. These signals will be called frequency bins (holding the energy of their respective frequency bands) and indexed with k. The coefficients inside each bin will be indexed by m. After the creation of the packet table, a basis is selected that emphasises the differences between the classes. This difference is measured with the help of a time-frequency energy map, defined as follows: Γc(j, k, m) = Nc X i (ˆxc i(j, k, m))2/ Nc X i ||xc i||2(1) where (j, k, m)denotes the position in the packet table, at splitting level j, frequency band kand coefficient mwithin the bin; ˆxc i(j, k, m)denotes the wavelet coefficient of click sample iand class cat position (j, k, m);xc ithe click sample iof class c;Ncthe number of training samples in class c. This map basically sums the packet tables of the clicks within one class and allows the comparison of the energy in a specific bin (j, k, ·)between different classes. This discrepancy can be measured with the following function, D(j, k, ·) = X m C−1 X p=1 C X q=p+1 D(Γp(j, k, m),Γq(j, k, m)) Here the difference between every pair of classes is measured through an additive discriminant function D, for which we used the squared l2-norm. A high value for Dmeans that that specific bin may be able to separate at least 2 classes that lie far apart. The local basis can now be selected using the following rule, if the measure on a bin D(j, k, ·), is higher than the sum of the measures over the two bins it splits into, D(j+ 1,2k, ·) + D(j+ 1,2k+ 1,·), then it is selected, otherwise it is split. After the local discriminant basis was created, we selected the 15 strongest coefficients, according to Fisher’s Algorithms 2009,21236 discriminant given by : FD =Pc(sc i−meanc(sc i))2 Pcvari(sc i)(2) where sare coefficients taken from a specific entry in the discriminating basis, and both the bar and varitake the mean and variance over all samples siin class cand meanctakes the mean over all classes. Essentially, this expression measures the distance between the class means and their common centre with respect to their widths, leading to high values when samples in a class lie tightly around their class centre. The 15 strongest features found using the selection procedure described above are summarised in Table 1. It is apparent that the best features were found in the lowest frequencies, where a dominant frequency can be found. Another reason why this may happen is that higher frequencies seem to be more variable in time, as was shown in Figure 1. A shift of a high frequency component will also shift its corresponding energy in the wavelet filter bank. The discrete wavelet transform with a symmlet wavelet as was used for this report has a very slow roll-off. After the down sampling step, aliasing will appear in and propagate through the filter bank. Normally, the aliasing (and phase distortion) is repaired in the synthesis step. This behaviour is not a problem for classification as long as it is consistent. However, phase delays in higher frequencies for consecutive clicks that cause subsequent energy shifts in the wavelet coefficients will act as a source of noise. Therefore, after it was confirmed that there were no high frequency features of interest, the algorithm was focussed on the lowest frequencies using a fifth order Butterworth low-pass filter at 2,000 Hz. After filtering the clicks were down sampled to remove redundant information. When this classification approach is followed on other data sets, the discriminant features should probably always first be selected from the full bandwidth to ensure optimal classification. Table 1. Selected wavelet packet coefficients for discrimination. The given index (k, m) is the position mof the coefficient in frequency band kat the given split level. The splitting level starts at 1, which indicates the original signal, i.e. at level 5 the signal was filtered 4 times. index split frq band (Hz) power index split frq band (Hz) power (1,2) 5 1 - 1,500 3.5 (1,1) 5 1 - 1,500 3.4 (1,3) 5 1 - 1,500 3.0 (1,4) 5 1 - 1,500 2.3 (2,11) 5 1,500 - 3,000 1.3 (1,31) 5 1 - 1,500 1.2 (1,9) 5 1 - 1,500 1.0 (1,21) 5 1 - 1,500 0.95 (1,10) 5 1 - 1,500 0.80 (2,8) 5 1,500 - 3,000 0.80 (1,13) 5 1 - 1,500 0.79 (1,18) 5 1 - 1,500 0.79 (1,14) 5 1 - 1,500 0.76 (1,22) 5 1 - 1,500 0.73 (1,16) 5 1 - 1,500 0.65 - All classifiers were trained using the characteristics of the first 50 clicks of the seven sets, reflecting the situation where the start of a recording is manually separated by an expert (or automatically when Algorithms 2009,21237 possible), and the remaining data would be processed automatically by the computer. The other clicks within the sets, and the eighth set, were then used for validation. Figure 2. Scatter plots of the four most discriminating features for five animals. The combination of these four characteristics already shows possible separation of five animals. Moreover, the features show a strong clustering tendency that suggests the use of a Gaussian model. Feature 1 Feature 2 Feature 3 Feature 4 Figure 3. Variability in the two strongest features (first feature on top, second on bottom) from Table 1during the dive. 0 200 400 600 800 1000 1200 −0.5 0 0.5 Click number Magnitude 0 200 400 600 800 1000 1200 −0.5 0 0.5 Click number Magnitude Figure 2shows a combination of two graphics with the four strongest features, measured using Fisher’s power of discrimination with Equation (2), from five animals that allows some insight in the feature space. The combination of just these four features already shows some possibility of separating the animals. Three animals are already separated in the left figure, while two are completely mixed. Combination with the two other characteristics in the right figure allows these two to be separated as well. An important observation in these figures, and one that we used to design the classifier, was that the data showed a strong clustering tendency. This suggested the application of a RBF network which has a natural way of modelling these clusters in its hidden layer, or a SVM which does not model the clusters themselves but their borders. Another reason we decided to use RBF and SVM based classifiers is because these methods have a local response defined by the distance of a sample to the cluster centres or support vectors. This is different from, for example, multi-layer perceptron networks, where a node in the first layer will give the same response for all points on a specific hyperplane. Although, as a result, perceptron networks may have better generalisation in areas of the feature space that are poorly sampled, Algorithms 2009,21238 we preferred the local properties that allowed more accurate modelling of the feature space as it is presented in Figure 2. To have an idea about the variability of the features, Figure 3shows the two strongest features during the whole dive. From this figure it can be suspected that 50 consecutive samples may not always be sufficient to characterise the variance of the feature. 3. Classification Description 3.1. Radial basis function network A detailed description about RBF networks can be found in [21]; its schematic is shown in Figure 4. Figure 4. Schematic of an RBF-network. An M-dimensional sample senters on the left, and is first run through the nhidden layer nodes where the distances between sand centres µare evaluated through Gaussian functions kn(s, µ). The outputs of the nGaussian functions are then weighed with weights wij and linearly combined in the second layer nodes (containing one node per class). Bias is usually represented with the help of a hidden layer node with a constant activation function, k0≡1. Each output layer then has a corresponding additional weight (w0j) that accounts for the bias factor. Taking the output vector y, the class of the sample sis computed by arg maxiyi. Input ' & $ % s1 . . . sM First layer º ¹· ¸ . . . . . . º ¹· ¸ ¶ µ³ ´ k1(s, µ) kn(s, µ) bias Second layer º ¹· ¸ . . . . . . º ¹· ¸ Piwi1ki(s, µ) Piwic ki(s, µ) Output ' & $ % y1 . . . yc An important reason to consider this type of network was that its two layers can be trained separately, without the use of a non-linear optimisation routine. This allows the training stage to be executed fast and makes it suitable for real-time applications. For the hidden layer activation functions we chose the Gaussian function given by k(s, µ) = exp(−||s−µ||2 2σ2)(3) where sis the input feature vector, µcontrols the function’s centre, and σits width. This layer was trained using a clustering algorithm, placing centres on top of dense locations in the data. One clustering algorithm that is used often is k-means [22], but this has the drawback that the number of clusters khas to be given in advance. There are various clustering methods that can search for an optimal (according to some defined statistic) number of clusters. One such method is described in [23], and we used a slightly adapted version. Initially, the clustering process starts with two clusters, splitting the data with Algorithms 2009,21239 k-means. A cluster was accepted and removed from the feature space when its projection in the direction of its principal component resembled a normal distribution, where normality was measured using the Anderson-Darling statistic [24,25]. The value for kwas then adjusted for the removed clusters and increased by one. The remaining clusters were combined and clustered again. This process was repeated until all clusters were accepted. To prevent random outcomes, k-means was always initialised by placing an additional centre in a cluster on the feature vector furthest away from the data’s centre [26]. In order to take advantage of class information in the clustering process, the clustering was done on the individual classes, instead of on all the data as a whole. Once the first layer has been defined by the clusters, training of the second layer is trivial; the number of nodes was set to the number of classes, using binary encoding for the targets (e.g., a sample from class 1 has target [1 0 0]t, and a sample from class 3 has target [0 0 1]t). Calculation of the weights, using a sum-of-squares error function, is then a fast linear process [21]. 3.2. Support vector machine classification Another popular network model that uses Gaussian functions are support vector machines (SVM). SVM generally solve the two class problem with a network structure that is similar to the RBF structure in Figure 4. Since there are only two classes, SVM only have one output. The main difference between the two networks relies on the underlying approach and in the way they are trained. Where the RBF network places Gaussian kernels on the cluster centres formed by the features in the feature space, the SVM network places the Gaussian kernels on those samples that define the boundary between two classes. These points are found by creating a separating hyperplane inside a higher dimensional feature space. A reason to consider SVM is that they tend to show strong generalisation performance [22]. In the case of SVM, the interpretation of the first layer is that it projects the data to a feature space where the two classes can be separated by a hyperplane. The position of the hyperplane is decided by those points that lie within a certain margin between the two classes (the support vectors). The projected feature space is allowed (or preferred) to have a higher dimension, as there is no penalty in the form of computational drawbacks or the ’curse of dimensionality’ since the actual mapping is never performed. All necessary calculations in the projected feature space are evaluated through inner products, which can be done with a kernel function. In our case this was the Gaussian function (3) that was also used for the RBF network. The training stage for SVM both defines the number of nodes in the hidden layer and calculates the output layer weights at the same time. There several approaches to train the network in the case of non-separable classes. First we looked at C-SVM, which aims at solving the following optimisation problem [27]: min w,ξ,b ||w||2+C N X i=1 ξi(4) under the constraint that the training data should classify correctly, and where wis the direction of the normal vector on the separating plane, Nthe number of samples, ξare the slack vectors which are nonzero for the samples that lie within the margin, and Cis the pre-defined penalty on these points. This problem has a straightforward optimal solution wousing Lagrange multipliers [27], expressed as follows Algorithms 2009,21240 wo= N X i=1 αitiφ(si)(5) where tiare the targets of mapped inputs si. Only a few of the Lagrange multipliers αiwill be non-zero, and those define both the first layer support vectors (si) and the second layer weights. It should be noted that wois never actually evaluated as the classification function given by f(s) = wT oφ(s) + bo=X i∈sv αitiφ(si)·φ(s) + bo(6) evaluates the inner product directly though the kernel function. A second training method that trains a SVM network regulating classification errors is ν-SVM [15]. The optimisation problem that is solved is then given by min w,ξ,b,ρ 1 2||w||2−νρ +1 N N X i=1 ξi(7) under the constraint of correct classification, and where ρis the margin in the mapped feature space and is left as a free variable. The fixed constant νgives some control over both the error rate on the training set (Pe, the ratio between misclassifications and number of training patterns) and the number of support vectors (Ns) from the optimisation through the relationships [22] : Pe≤ν;and Nν ≤Ns.(8) This can be a useful property as the number of support vectors play a large role on the speed of the network. The SVM network only separates two classes, but the same algorithm can easily be extended to several classes, for example by using a one-against-one approach. For this classification method a separate machine is created for every combination of two classes, meaning that when there are nclasses this results in n(n−1) 2machines. A new pattern is then presented to all the machines, and the class that occurs most frequently in the outcome is chosen. In the case when there are two or more classes with identical frequencies, the pattern is defined as unclassifiable. 4. Classification Results In order to classify the patterns, the 15 strongest features according to equation (2) were selected from the local discriminant basis. The C-SVM results were obtained using a toolbox from [28]; custom code in Matlab was written to obtain the ν-SVM results. For both SVM approaches a total number of 21 support vector machines was used to classify the seven classes; every support vector machine was trained with 100 patterns. Each machine has two parameters that can be tuned, the width and either the parameter Cor ν. To simplify the training stage, their values were kept identical over all machines. In order to select suitable parameters, the machines were trained while varying one parameter and keeping the other one fixed. At each value the classifier was trained 10 times with noise added to the data. The noise was drawn from a zero-mean normal distribution with the standard deviation taken from the standard deviation of the features. Algorithms 2009,21247 29. Halkias, X.; Ellis, D. Estimating the number of marine mammals using recordings from one microphone. In Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing ICASSP-06, Toulouse, France, May 14-19, 2006. 30. Ioup, J.; Ioup, G. Self-organizing maps for sperm whale identification. In Proceedings ot TwentyThird Gulf of Mexico Information Transfer Meeting; McKay, M., Nides, J., Eds.; U.S. Department of the Interior, Minerals Management Service: New Orleans, LA, USA, January 11, 2005; pp. 121-129. c °2009 by the authors; licensee Molecular Diversity Preservation International, Basel, Switzerland. This article is an open-access article distributed under the terms and conditions of the Creative Commons Attribution license http://creativecommons.org/licenses/by/3.0/.