Automatic large-scale classification of bird sounds is strongly improved by unsupervised feature learning

Spectral features and feature learning

Raw audio data is not generally suitable input to a classification algorithm: even if the audio inputs were constrained to a fixed duration (so that the data dimensionality was constant), the dimensionality of an audio signal (considered as a vector) would be extremely large, and would not represent sound in such a way that perceptually similar sounds would generally be near neighbours in the vector space. Hence audio data is usually converted to a spectrogram-like representation before processing, i.e., the magnitudes of short-time Fourier transformed (STFT) frames of audio, often around 10 ms duration per frame. (Alternatives to STFT which have been considered for bird sound classification include linear prediction (Fox, 2008), wavelets (Selin, Turunen & Tanttu, 2007) and chirplets (Stowell & Plumbley, in press)). An STFT spectrum indicates the energy present across a linear range of frequencies. This linear range might not reflect the perceptual range of a listener, and/or the range of frequencies at which the signal carries information content, so it is common to transform the frequency axis to a more perceptual scale, such as the Mel scale originally intended to represent the approximately logarithmic frequency sensitivity of human hearing. This also reduces the dimensionality of the spectrum, but even the Mel spectrum has traditionally been considered rather high-dimensional for automatic analysis. A convention, originating from speech processing, is to transform the Mel spectrum using a cepstral analysis and then to keep the lower coefficients (e.g., the first 13) which typically contain most of the energy (Davis & Mermelstein, 1980). These coefficients, the Mel frequency cepstral coefficients (MFCCs), became widespread in applications of machine learning to audio, including bird vocalisations (Stowell & Plumbley, 2010).

MFCCs have some advantages, including that the feature values are approximately decorrelated from each other, and they give a substantially dimension-reduced summary of spectral data. Dimension reduction is advantageous for manual inspection of data, and also for use in systems that cannot cope with high-dimensional data. However, as we will see, modern classification algorithms can cope very well with high-dimensional data, and dimension reduction always reduces the amount of information that can be made available to later processing, risking discarding information that a classifier could have used. Further, there is little reason to suspect that MFCCs should capture information optimal for bird species identification: they were designed to represent human speech, yet humans and birds differ in their use of the spectrum both perceptually and for production. MFCCs aside, one could use raw (Mel-)spectra as input to a classifier, or one could design a new transformation of the spectral data that would tailor the representation to the subject matter. Rather than designing a new representation manually, we consider automatic feature learning.

The topic of feature learning (or representation learning, dictionary learning) has been considered from many perspectives within the realm of statistical signal processing (Bengio, Courville & Vincent, 2013; Jafari & Plumbley, 2011; Coates & Ng, 2012; Dieleman & Schrauwen, 2013) . The general aim is for an algorithm to learn some transformation that, when applied to data, improves performance on tasks such as sparse coding, signal compression or classification. This procedure may be performed in a “supervised” manner, meaning it is supplied with data as well as some side information about the downstream task (e.g., class labels), or “unsupervised”, operating on a dataset but with no information about the downstream task. A simple example that can be considered to be unsupervised feature learning is principal components analysis (PCA): applied to a dataset, PCA chooses a linear projection which ensures that the dimensions of the transformed data are decorrelated (Bengio, Courville & Vincent, 2013). It therefore creates a new feature set, without reference to any particular downstream use of the features, simply operating on the basis of qualities inherent in the data.

Recent work in machine learning has shown that unsupervised feature learning can lead to representations that perform very strongly in classification tasks, despite their ignorance of training data labels that may be available (Coates & Ng, 2012; Bengio, Courville & Vincent, 2013). This rather surprising outcome suggests that feature learning methods emphasise patterns in the data that turn out to have semantic relevance, patterns that are not already made explicit in the basic feature processing such as STFT. A second surprising aspect is that such representations often perform the opposite of feature reduction, increasing the dimensionality of the problem without adding any new information: a deterministic transformation from one feature space to a higher-dimensional feature space cannot, in an information-theoretic sense, add any information that is not present in the original space. However, such a transformation can help to reveal the manifold structure that may be present in the data (Olshausen & Field, 2004) . Neural networks, both in machine implemetations and in animals, perform such a dimension expansion in cases where one layer of neurons is connected as input to a larger layer of neurons (Olshausen & Field, 2004).

In our study, however, we will not use a feature learning procedure intended to parallel a biological process. Instead, we use spherical k-means, a simple and highly scalable modification of the classic k-means algorithm (Coates & Ng, 2012; Dieleman & Schrauwen, 2013). We perform a further adaptation of the algorithm to ensure that it can run in streaming fashion across large audio datasets, to be described in ‘Materials and Methods’.

Birdsong often contains rapid temporal modulations, and this information should be useful for identifying species-specific characteristics (Stowell & Plumbley, in press). From this perspective, a useful aspect of feature learning is that it can be applied not only to single spectral frames, but to short sequences (or “patches”) of a few frames. The representation can then reflect not only characteristics of instantaneous frequency patterns in the input data, but characteristics of frequencies and their short-term modulations, such as chirps sweeping upwards or downwards. This bears some analogy with the “delta-MFCC” features sometimes used by taking the first difference in the time series of MFCCs, but is more flexible since it can represent amplitude modulations, frequency modulations, and correlated modulations of both sorts (cf. Stowell & Plumbley, in press). In our study we tested variants of feature learning with different temporal structures: either considering one frame at a time (which does not capture modulation), multiple frames at a time, or a variant with two layers of feature learning, which captures modulation across two timescales.

Datasets

We gathered four datasets, each representing a large amount of audio data and a large number of species to classify (Table 1). Two of the datasets (nips4b and lifeclef2014) consist of the publicly-released training data from bird classification challenges organised by the SABIOD project (Glotin et al., 2013; Goëau et al., 2014). The nips4b dataset is multilabel (median 1 species per recording, range 0–6); the lifeclef2014 dataset is single-label but much larger. (Some of the data in lifeclef2014 includes annotations of “background species” which could be used alongside the primary annotation to construct a multilabel task; we did not do this.) Note that we only use the publicly-released training data from those challenges, and not any private test data, and so our evaluation will be similar in nature to their final results but not precisely comparable. For evaluation we partitioned each of these datasets into two, so that we could run two-fold crossvalidation: training on one half of the dataset and testing on the other half, and vice versa.

In addition, the British Library Sound Archive has a large collection of environmental sound recordings, and they made available to us a subset of 60 “dawn chorus” recordings. This consisted of 20 recordings each from three UK-based recordists, ranging in duration from 2 min to 20 min, and annotated by each recordist with a list of species heard (median 6 species per recording, range 3–12). We refer to this dataset as bldawn, and perform three-fold stratified crossvalidation: for each recordist, we train the system using the data from the other two recordists, and then test on the audio from the held-out recordist. This stratified approach is useful because it tests whether the system can generalise to recordings from unknown recordists, rather than adapting to any specifics of the known recordists.

We also gathered a single-label dataset as a subset of the recordings available from the Xeno Canto website¹, covering many of the common UK bird species, and covering at least all the species present in the bldawn dataset. We refer to this dataset as xccoverbl. For each species included, we queried Xeno Canto to retrieve three different recordings, preferring to retrieve recordings from the UK, but allowing the system to return recordings from further afield if too few UK recordings were available. Our search query also requested high-quality recordings (quality label ‘A’), and song rather than calls, where possible. Since we retrieved three examples for each species, this enabled us to partition the dataset for three-fold crossvalidation: not stratified into individual recordists (as was bldawn), but sampled from a wide range of recordists.

These datasets have widely varying characteristics, for example in the typical duration of the sound files, the recording location, and the number of classes to distinguish (Table 1). Note that most of the datasets have different and irreconcilable lists of class labels: in particular, for bldawn and xccoverbl the class label is the species, whereas nips4b and lifeclef2014 use separate labels for song and calls. Of our datasets only bldawn and xccoverbl have strong overlap in their species lists. Therefore only these datasets could be combined to create larger pools of training data.

In this work we performed automatic classification for each audio file, without any segmentation procedure to select region(s) of bird vocalisation in the file. The only segmentation that is done is implicit in the collection processes for the dataset: for the two datasets originating from Xeno Canto, each audio clip might or might not contain a large amount of silence or other noise, depending on the contributor; for nips4b the audio is collected from remote monitoring stations with no manual selection; for bldawn the audio is selected by the contributor, but not trimmed to a specific vocalisation, instead selected to present a long dawn chorus audio recording.

Feature learning method

As discussed in ‘spectral features and feature learning’, the aim of unsupervised feature learning is to find some transformation of a dataset, driven only by the characteristics inherent in that dataset. For this we use a method that has shown promise in previous studies, and can be run effectively at big data scales: spherical k-means, described by Coates & Ng (2012) and first applied to audio by Dieleman & Schrauwen (2013). There are many feature-learning methods available, including neural networks such as restricted Boltzmann machines (used e.g., in Erhan et al., 2010), or methods based on sparse coding such as K-SVD (Aharon, Elad & Bruckstein, 2006). Our choice of method is motivated by the promising results of Dieleman & Schrauwen (2013) but also by our imperative to enable feature learning at very large scale. This leads to a preference for techniques of low computational complexity, and which can be applied to data in streaming fashion.

Spherical k-means is related to the simple and well-known k-means clustering algorithm (Lloyd, 1982), except that instead of searching for cluster centroids which minimise the Euclidean distance to the data points, we search for unit vectors (directions) to minimise their angular distance from the data points. This is achieved by modifying the iterative update procedure for the k-means algorithm: for an input data point, rather than finding the nearest centroid by Euclidean distance and then moving the centroid towards that data point, the nearest centroid is found by cosine distance, (1)${\displaystyle \text{cosine distance}=1-cos\left(\theta \right)=1-\frac{A\cdot B}{\Vert A\Vert \Vert B\Vert }\hspace{0.17em},}$ where A and B are vectors to be compared, θ is the angle between them, and ‖⋅‖ is the Euclidean vector norm. The centroid is renormalised after update so that it is always a unit vector. Figure 1 shows an example of spherical k-means applied to synthetic data. Spherical k-means thus finds a set of unit vectors which represent the distribution of directions found in the data: it finds a basis (here an overcomplete basis) so that data points can in general be well approximated as a scalar multiple of one of the basis vectors. This basis can then be used to represent input data in a new feature space which reflects the discovered regularities, in the simplest case by representing every input datum by its dot product with each of the basis vectors (Coates & Ng, 2012; Dieleman & Schrauwen, 2013): (2)${\displaystyle x^{\prime }\left(n,j\right)=\sum _{i=1}^M{b}_j\left(i\right)x\left(n,i\right),}$ where x represents the input data indexed by time frame n and feature index i (with M the number of input features, e.g., the number of spectral bins), b_j is one of the learnt basis vectors (indexed by j∈[1, k]), and x′ is the new feature representation. In our case, the data on which we applied the spherical k-means procedure consisted of Mel spectral frames (M = 40 dimensions), which we first normalised and PCA-whitened as in Dieleman & Schrauwen (2013).

We also tested configurations in which the input data was not one spectral frame but a sequence of them—e.g., a sequence of four spectral frames at a time—allowing the clustering to respond to short-term temporal patterns as well as spectral patterns. We can write this as (3)${\displaystyle x^{\prime }\left(n,j\right)=\sum _{\delta =0}^{\Delta -1}\sum _{i=1}^M{b}_j\left(\delta ,i\right)x\left(n+\delta ,i\right)\hspace{0.17em},}$ where Δ is the number of frames considered at a time, and the b are now indexed by a frame-offset as well as the feature index. (See Fig. 10 to preview examples of such bases.) Alternatively, this can be thought of as “stacking” frames, e.g., stacking each sequence of four 40-dimensional spectral frames to give a 160-dimensional vector, before applying (2) as before. In all our experiments we used a fixed k = 500, a value which has been found useful in previous studies (Dieleman & Schrauwen, 2013).

The standard implementation of k-means clustering requires an iterative batch process which considers all data points in every step. This is not feasible for high data volumes. Some authors use “minibatch” updates, i.e., subsamples of the dataset. For scalability as well as for the potential to handle real-time streaming data, we instead adapted an online streaming k-means algorithm, “online Hartigan k-means” (McFee, 2012, Appendix B). This method takes one data point at a time, and applies a weighted update to a selected centroid dependent on the amount of updates that the centroid has received so far. We adapted the method of (McFee, 2012, Appendix B) for the case of spherical k-means. k-means is a local optimisation algorithm rather than global, and may be sensitive to the order of presentation of data. Therefore in order to minimise the effect of order of presentation for the experiments conducted here, we did not perform the learning in true single-pass streaming mode. Instead, we performed learning in two passes: a first streamed pass in which data points were randomly subsampled (using reservoir sampling) and then shuffled before applying PCA whitening and starting the k-means procedure, and then a second streamed pass in which k-means was further trained by exposing it to all data points. Our Python code implementation of online streaming spherical k-means is available as Supplemental Information.

As a further extension to the method, we also tested a two-layer version of our feature-learning method, intended to reflect detail across multiple temporal scales. In this variant, we applied spherical k-means feature learning to a dataset, and then projected the dataset into that learnt space. We then downsampled this projected data by a factor of 8 on the temporal scale (by max-pooling, i.e., taking the max across each series of 8 frames), and applied spherical k-means a second time. The downsampling operation means that the second layer has the potential to learn regularities that emerge across a slightly longer temporal scale. The two-layer process overall has analogies to deep learning techniques, most often considered in the context of artificial neural networks (Erhan et al., 2010; Bengio, Courville & Vincent, 2013), and to the progressive abstraction believed to occur towards the higher stages of auditory neural pathways.

Classification and evaluation

Our full classification workflow started by converting each audio file to a standard sample-rate of 44.1 kHz. We then calculated Mel spectrograms for each file, using a frame size of 1024 frames with Hamming windowing and no overlap. We chose no overlap rather than 50% overlapped frames simply to reduce the volume of data to be processed. We filtered out spectral energy below 500 Hz, a heuristic choice which strongly reduces the amount of environmental noise present, a benefit which is traded off against the cost that this will discard some energy from species that vocalise below 500 Hz (such as the eagle-owl Bubo bubo). We then normalised the root-mean-square (RMS) energy in each spectrogram.

For each spectrogram we then optionally applied the noise-reduction procedure that we had found to be useful in our NIPS4B contest submission (Stowell & Plumbley, 2013), a simple and common median-based thresholding. This consists of finding the median value for each spectral band in a spectrogram, then subtracting this median spectrum from every frame, and setting any resulting negative values to zero. This therefore preserves only the spectral energy that rises above the median bandwise energy. In principle it is a good way to reduce the stationary noise background (Fig. 2), but is not designed to cope well with fluctuating noise. However its simplicity makes it easy to apply across large datasets efficiently.

The Mel spectrograms, either noise-reduced or otherwise, could be used directly as features. We also tested their reduction to MFCCs (including delta features, making 26-dimensional data), and their projection onto learned features, using the spherical k-means method described above. For the latter option, we tested projections based on single frame as well as on sequences of 2, 3, 4 and 8 frames, to explore the benefit of modelling short-term temporal variation. We also tested the two-layer version based on the repeated application to 4-frame sequences across two timescales.

The feature representations thus derived were all time series. In order to reduce them to summary features for use in the classifier, we tested two common and simple techniques: summarising each feature dimension independently by its mean and standard deviation, or alternatively by its maximum. These are widespread but are not designed to reflect any temporal structure in the features (beyond the fine-scale temporal information that is captured by some of our features). Therefore, for the Mel and MFCC features we also tested summarising by modulation coefficients: we took the short-time Fourier transform (STFT) along the time axis of our features, and then downsampled the spectrum to a size of 10 to give a compact representation of the temporal evolution of the features (cf. Lee, Han & Chuang, 2008). The multi-frame feature representations already intrinstically included short-term summarisation of temporal variation, so to limit the overall size of the experiment, for the learned feature representations we only applied the mean-and-standard-deviation summarisation. Overall we tested six types of non-learned representation against six types of learned representation (Table 2).

To perform classification on our temporally-pooled feature data, then, we used a random forest classifier (Breiman, 2001). A random forest classifier is an ensemble method which trains many decision-tree classifiers on the same dataset: the decision trees are different from each other due to the use of “bagging”—drawing a different bootstrap sample from the training dataset for each tree—and also by considering only a small random subset of the available data features as candidates for each split. This randomisation reduces the correlation between individual decision trees. To make a prediction, the random forest uses a simple vote to aggregate the predictions of its decision trees: in this work we use probabilistic outputs from the classifier, meaning that the vote proportions are reinterpreted as probabilities. Random forests and other tree-ensemble classifiers perform very strongly in a wide range of empirical evaluations (Caruana & Niculescu-Mizil, 2006), and were used by many of the strongest-performing entries to the SABIOD evaluation contests (Glotin et al., 2013; Fodor, 2013; Potamitis, 2014). For this experiment we used the implementation from the Python scikit-learn project (Pedregosa et al., 2011). Note that scikit-learn v0.14 was found to have a specific issue preventing training on large data, so we used a pre-release v0.15 after verifying that it led to the same results with our smaller datasets.

We did not manually tune any parameters of the classifier: parameter tuning can lead to improvements in performance, but can also lead to overfitting to particular dataset characteristics, so in all cases we trained a random forest with 200 trees using the ‘entropy’ (information-gain) criterion to measure the quality of a split. However, since our experiment covered both single-label and multilabel classification, we did test three different ways of using the classifier to make decisions:

1. Single-label classification: this assumes that there is only one species present in a recording. It therefore cannot be applied to multilabel datasets, but for single-label datasets it may benefit from being well-matched to the task.

2. Binary relevance: this divides the multilabel classification task into many single-label tasks, training one separate classifier for each of the potential output labels (Tsoumakas, Katakis & Vlahavas, 2010). This strategy ignores potential correlations between label occurrence, but potentially allows a difficult task to be approximated as the combination of more manageable tasks. Binary relevance is used e.g., by Fodor (2013).

3. Full multilabel classification: in this approach, a single classifier (here, a single random forest) is trained to make predictions for the full multi-label situation. Predicting presence/absence of every label simultaneously can be computationally difficult compared against a single-label task, and may require larger training data volumes, but represents the full situation in one model (Tsoumakas, Katakis & Vlahavas, 2010).

For each of these methods the outputs from the classifier are per-species probabilities. We tested all of our datasets using the full multilabel classifier, then for comparison we tested the single-label datasets using the single-label classifier, and the multi-label datasets using the binary-relevance classifier.

Some of our datasets contain long audio recordings, yet none of the annotations indicate which point(s) in time each species is heard. This is a common format for annotations: for example, the bldawn annotations are derived directly from the archival metadata, which was not designed specifically for automatic classification. Long audio files present an opportunity to make decisions either for the entire file as one scene, or in smaller “decision windows”, for which the decisions are then pooled to yield overall decisions. We tested this empirically, using decision windows of length 1, 5 or 60 s or the whole audio. Each decision window was treated as a separate datum for the purposes of training and testing the classifier, and then the decisions were aggregated per audio file using either mean or maximum. The mean probability of a species across all the decision windows is a reasonable default combination; we compared this against the maximum with the motivation that if a bird is heard only at one point in the audio, and this leads to a strong detection in one particular decision window, then such a strong detection should be the overriding factor in the overall decision. For some datasets (nips4b) we did not test long windows since all audio files were short; while for lifeclef2014 we used only whole-audio classification because of the runtime costs of evaluating these combinations over this largest dataset.

We performed feature learning, training and testing separately for each of our four datasets, using the appropriate two- or threefold crossvalidation described above, and across all combinations of the feature settings we have just described. Figure 3 summarises the main stages of the workflow described.

We evaluated the performance in each experimental run using two measures, the Area Under the ROC Curve (AUC) and the Mean Average Precision (MAP). The AUC statistic is an evaluation measure for classification/detection systems which has many desirable properties (Fawcett, 2006): unlike raw accuracy, it is not affected by “unbalanced” datasets having an uneven mixture of true-positive and true-negative examples; and it has a standard probabilistic interpretation, in that the AUC statistic tells us the probability that the algorithm will rank a random positive instance higher than a random negative instance. Chance performance is always 50% for the AUC statistic. The AUC is a good statistic to use to evaluate the output from a probabilistic classifier in general, but for user-facing applications, which may for example show a ranked list of possible hits, the Mean Average Precision (MAP) statistic leads to an evaluation which relates more closely to user satisfaction with the ranked list. The key difference in effect is that the AUC statistic applies an equal penalty for a misordering at any position on the ranked list, whereas the MAP statistic assigns greater penalties higher up the ranking (Yue et al., 2007). We therefore calculate both evaluation statistics.

To test for significant differences in the performance statistics, we applied a generalised linear model (GLM) using the lme4 package (Bates et al., 2014) for R 2.15.2 (R Core Team, 2012). We focused primarily on the AUC for significance testing, since the AUC and MAP statistics are related analyses of the same data. Since AUC is bounded in the range [0, 1], we applied the GLM in the logistic domain: note that given the probabilistic interpretation of AUC, the logistic model is equivalent to applying a linear GLM to odds ratios, a fact which facilitates interpretation. Every experimental run for a given dataset used the same set of folds, so we used a repeated-measures version of the GLM with the “fold index” as the grouping variable. We tested for individual and pairwise interactions of our five independent categorical variables, which were as follows:

• choice of feature set and temporal summarisation method, testing the 12 configurations listed in Table 2:

• noise reduction on vs. off;

• classifier mode (multilabel vs. either single-label or binary-relevance);

• decision pooling window duration (1, 5, or 60 s or whole audio);

• decision pooling max vs. mean.

Combinatorial testing of all these configurations resulted in 12 × 2 × 2 × 4 × 2 = 384 crossvalidated classification experiments for the bldawn and xccoverbl datasets. For the other datasets we did not test all four pooling durations, for reasons given above: the number of crossvalidated experiments was thus 192 for nips4b and 96 for lifeclef2014. Since the tests of lifeclef2014 did not vary decision pooling, decision pooling factors were not included in that GLM. We considered effects to be significant when the 95% confidence interval calculated from the GLM excluded zero, in which cases we report the estimated effects as differences in odds-ratios.

Additional tests

The bldawn dataset has relatively few annotations, since it only consists of 60 items. We therefore wanted to explore the use of auxiliary information from other sources to help improve recognition quality, in particular using the xccoverbl dataset, which has strong overlap in the list of species considered. In further tests we tested three ways of using this additional data:

1. Cross-condition training, meaning training on one dataset and testing on the other. The two datasets have systematic differences—for example, xccoverbl items are annotated with only one species each, and are generally shorter—and so we did not expect this to yield very strong results.

2. Data augmentation for the feature learning step, meaning that feature learning is conducted using the training data for the bldawn as well as all of the xccoverbl data. This gives a larger and more varied pool of data for the feature learning step, which we expected to give a slight improvement to the results of feature learning.

3. Data augmentation for feature learning and also for training. Although the systematic differences mean the xccoverbl training data might not guide the classifier in the correct way, it holds many more species annotations for the species of interest, in a wider set of conditions. We expected the combined training would provide stronger generalisation performance.

We evaluated these train/test conditions as separate evaluation runs.

We also wanted to distinguish between two possible explanations for any difference between the performance of the different feature sets: was it due to intrinsic characteristics of the features, or more simply due to the dramatic differences in feature dimensionality (which ranged 26–1,000; see Table 2)? Differences in dimensionality might potentially give more degrees of freedom to the classifier without necessarily capturing information in a useful way. In order to test this, we ran a version of our test in which for each experimental run we created a random projection matrix which projected the feature set to a fixed dimensionality of 200. For MFCC/Mel features this was a simple form of data expansion, while for learned features it was a form of data reduction. By standardising the feature dimensionality, this procedure decoupled the nature of the feature set from the degrees of freedom available to the classifier. We ran this test using the nips4b dataset.