# Load the necessary packages

Make sure you have the latest version of the packages:

    install.packages("bioacoustics")

    # The bioacoustics package may also be installed from GitHub using devtools as follows:
    install.packages("devtools")
    devtools::install_github("wavx/bioacoustics") # For the latest, unstable version

    install.packages("warbleR")
    install.packages("randomForest")

    # Load the packages
    library(warbleR)
    library(bioacoustics)
    library(tools)
    library(randomForest)

# Load audio files

We can use the `quer_xc()` function from warbleR to download bird vocalizations from Xeno-Canto: https://www.xeno-canto.org/

We are going to choose calls from Catharus bicknelli, songs from Passerella iliaca and Setophaga magnolia recorded in the United States and Canada.

We will filter only “A” quality recordings, then, pick up only the first nine, and merge all the metadata into a single “df” data frame. This data frame will be used to download MP3 files in your working directory directly from Xeno-Canto with the `quer_xc()` function:

    df1 = quer_xc(qword ='Catharus bicknelli type:call cnt:"United States"', download = FALSE)
    df1 = df1[df1$Vocalization_type=="call",]
    df1 = df1[df1$Quality=="A",]
    df1 = df1[1:9,]

    df2 = quer_xc(qword ='Setophaga magnolia type:song cnt:"Canada"', download = FALSE)
    df2 = df2[df2$Quality=="A",]
    df2 = df2[1:9,]

    df3 = quer_xc(qword ='Passerella iliaca type:song cnt:"Canada"', download = FALSE)
    df3 = df3[df3$Vocalization_type=="song",]
    df3 = df3[df3$Quality %in% c("A", "B"),]
    df3 = df3[1:9,]

    df = rbind(df1,df2,df3)
    rm(df1,df2,df3)

    # Visualize your data frame
    View(df)

    # Create data directory if not exists
    if(!dir.exists("data"))
      dir.create("data")

    # Download the MP3 files into your data directory
    quer_xc(X = df, download = TRUE, path = "data")

Now that we have recordings stored in the data directory, we can read one of them to look at its structure and content. Let's use the `read_audio()` function in bioacoustics to manually read a recording of Catharus bicknelli:

    CATBIC <- read_audio("data/Catharus-bicknelli-54864.mp3")
    CATBIC

We can see that the MP3 file has been converted into a Wave object with 7551360 samples, a duration of 157.32 seconds, a sampling rate of 48000 Hz, a bit depth of 16 bits, and that it contains one channel (mono, stereo being two channels).

Remember that you just have to divide the number of samples by the sampling rate to retrieve the duration (s) of an audio file.

We will use the Fast Fourrier Transform (FFT) to compute a frequency-time representation of the recording called spectrogram. A spectrogram being the representation of the spectrum of frequencies in a recording as they vary through time. This representation is used to detect animal vocalizations and extract acoustic features useful for classification, with the purpose of animal identification.

# Plot audio files

There are several options to display animal vocalizations in audio files with R. You can use both `spectro()` or `fspec()` functions to generate spectrograms with bioacoustics. `fspec()` generates only a matrix of the spectrogram, and thus has to be used with the `image()` function to display the spectrogram. It is also possible to use the `spectro()` function in Seewave.

Next, we will search manually and display an audio event (here, a bird vocalization) from a recording of Catharus bicknelli. To display a Region Of Interest (ROI) of the recording we will use temporal and frequency filters. Let's start with a temporal slice from 1 to 10 secs and a FFT size of 512 samples.

    # Set plot margins to 0
    par(mar = c(0, 0, 0, 0), oma = c(0, 0, 0, 0))

    # Display with spectro()
    ticks <- c(from = 1000, to = 20000, by = 1000) # frequency tick marks from 1 to
                                                   # 20 kHz, and steps at each 1 kHz
    temp_slice <- c(from = 1, to = 10) # in seconds
    spectro(CATBIC, tlim = temp_slice, FFT_size = 512, ticks_y = ticks)

Let's display spectrograms with various time / frequency limits (with `tlim=` and `flim=` arguments). You can also play with other arguments in `spectro()` and `fspec()` functions such as the percent of overlap between two FFTs (with `FFT_overlap=`) and various FFT resolutions (with `FFT_size=`). Note that the arguments are briefly explained in the documentation of each function:

    # Access the arguments of the spectro function
    ?spectro
    ?fspec

First, let's shorten the temporal axis from 2 to 3.5 secs to work on a shorter time window and compare the outputs from `spectro()` and `fspec()` functions. Note that spectrograms can also be generated automatically while using the detection functions from bioacoustics. We will explore that in details in section 4.1.

    # Set plot margins to 0
    par(mar = c(0, 0, 0, 0), oma = c(0, 0, 0, 0))

    # Display the spectrogram with spectro()
    ticks <- c(from = 1000, to = 20000, by = 1000) # frequency tick marks from 1 to
                                                   # 20 kHz, 1 kHz steps
    temp_slice <- c(from = 2, to = 3.5) # in seconds
    spectro(CATBIC, tlim = temp_slice, FFT_size = 512, ticks_y = ticks)

    # fspec() gives you the spectrogram matrix with energy values (dB)
    spec_mx <- fspec(CATBIC, tlim = temp_slice, FFT_size = 512, rotate = TRUE)

    # You can display the spectrogram with image()
    image(spec_mx, xaxt = "n", yaxt = "n")

The tick marks on the (frequency) y-axis were defined in the `spectro()` function starting from 1 to 20 kHz with an interval at each 1 kHz. The FFT size was 512 samples with an overlap between two FFT windows set by default at 0.875. Now try these settings: FTT size = 256, 1024 and 2048; FFT overlap = 0.3, 0.6, 0.9…

Another interesting thing to perform with the fspec outputs, is to implement your own set of filters. Let's try to reduce the background noise from a spectrogram with a narrower time and frequency bandwidth:

    temp_slice <- c(from = 2.5, to = 3.5)
    freq_slice <- c(from = 1500, to = 20000)
    spec_o <- fspec(CATBIC, tlim = temp_slice, flim = freq_slice, FFT_size = 512, rotate = TRUE)

    ## min and max (range) dB intensity
    range(spec_o) # -120 (min) to 0 dB (max)
    # Note that the tolerance of your recorders depends on the number of bits.
    # 16-bit recorders offer only around -96 dB tolerance and sound pressure above
    # this level is clipped to 0 dB.

    ## Let's try a filter by mean + sd intensity
    spec_f <- fspec(CATBIC, tlim = temp_slice, flim = freq_slice, FFT_size = 512, rotate = TRUE)
    spec_f[spec_f < mean(spec_f) + sd(spec_f)] <- -120
    # Works well with high intensity audio events, but leads to
    # false negatives (missed events) otherwise.

    par(mar = c(0, 0, 0, 0), oma = c(0, 0, 0, 0))
    image(spec_o, xaxt="n", yaxt="n")
    image(spec_f, xaxt="n", yaxt="n")

# The use of filters to detect and extract audio events

The functions used to detect and extract audio events in a recording also rely on “generic” filters based on frequency and duration, along with other “specific” filters. Let's take a quick look at the generic filters available in the `threshold_detection()` and `blob_detection()` functions:

- High Pass (HPF) and Low Pass filters (LPF) can be employed to reduce the amount of unwanted noise in the recording or to track particular audio events within a narrower frequency bandwidth than the recording sampling rate. Frequencies below the HPF and above the LPF cutoff are greatly attenuated. These frequency filters can be set using the `HPF=` and `LPF=` arguments in the `threshold_detection()` and `blob_detection()` functions.
- Minimum and maximum duration of an audio event, and a minimum time between two audio events also help reduce the amount of unwanted noise, or track a particular audio event within a narrower temporal window. These temporal filters can be set using the `min_dur=`, `max_dur=`, and `TBE=` arguments in both `threshold_detection()` and `blob_detection()` functions.
- Other set of filters are specific to each detection function and will be defined while working with these functions on bird vocalizations.

## Detect and extract audio events in a recording

### Threshold detection

Let's start with the `threshold_detection()` function on a recording containing calls from Catharus bicknelli. This function is an amplitude threshold detector that picks up audio events above the Signal to Noise Ratio (SNR). It combines several algorithms for detection, filtering and audio feature extraction. We will play with the arguments of this function to understand their implication in the detection and extraction of audio events (here, calls of Catharus bicknelli).

    # Access the arguments of the threshold_detection function
    ?threshold_detection

    # Set each argument according to the targeted audio events
    TD <- threshold_detection(
      CATBIC, # Either a path to an audio file (see ?read_audio), or a Wave object
      threshold = 12, # 12 dB SNR sensitivity for the detection algorithm
      time_exp = 1, # Time expansion factor of 1. Only needed for bat recordings.
      min_dur = 140, # Minimum duration threshold of 140 milliseconds (ms)
      max_dur = 440, # Maximum duration threshold of 440 ms
      TBE = 10, # Minimum time window between two audio events of 10 milliseconds
      EDG = 0.996, # Temporal masking with Exponential Decay Gain from 0 to 1
      LPF = 10000, # Low-Pass Filter of 10 kHz
      HPF = 1000, # High-Pass Filter of 1 kHz
      FFT_size = 256, # Size of the Fast Fourrier Transform (FFT) window
      FFT_overlap = 0.875, # Percentage of overlap between two FFT windows

      start_thr = 25, # 25 dB threshold at the start of the audio event
      end_thr = 30, # 30 dB threshold at the end of the audio event
      SNR_thr = 10, # 10 dB SNR threshold at which the extraction of the audio event stops
      angle_thr = 45, # 45° of angle at which the extraction of the audio event stops
      duration_thr = 440, # Noise estimation is resumed after 440 ms
      NWS = 1000, # Time window length of 1 s used for background noise estimation
      KPE = 1e-05, # Process Error parameter of the Kalman filter (for smoothing)
      KME = 1e-04, # Measurement Error parameter of the Kalman filter (for smoothing)

      settings = FALSE, #  Save on a list the above parameters set with this function
      acoustic_feat = TRUE, # Extracts the acoustic and signal quality parameters
      metadata = FALSE, # Extracts on a list the metadata embedded with the Wave file
      spectro_dir = "./Spectros", # Directory where to save the spectrograms
      time_scale = 1, # Time resolution of 2 ms for spectrogram display
      ticks = TRUE # Tick marks and their intervals are drawn on the y-axis (frequencies)
    )

    # Get the number of extracted audio events
    nrow(TD$event_data)

Let the HTML page open with the 56 spectrograms (each representing an extracted audio event). These settings will be our benchmark for the number of audio events that can be extracted with the `threshold_detection()` function. In the following exercise, you will try to reach or beat this number by exploring different combinations of parameters for each argument of the function.

    # Let's try various settings, starting with 1024 FFT size instead of 256.
    TD <- threshold_detection(
      CATBIC, threshold = 12, time_exp = 1, min_dur = 140, max_dur = 440, TBE = 10,
      EDG = 0.996, LPF = 10000, HPF = 1000, FFT_size = 1024, FFT_overlap = 0.875,
      start_thr = 25, end_thr = 30, SNR_thr = 10, angle_thr = 45, duration_thr = 440,
      NWS = 1000, KPE = 1e-05, KME = 1e-04, settings = FALSE, acoustic_feat = TRUE,
      metadata = FALSE, spectro_dir = "./Spectros", time_scale = 1,
      ticks = c(1000, 10000, 1000) # Tick marks from 1 to 10 kHz with 1 kHz interval
    )

    # Take a look at the spectrograms and compare them with the previous extraction.
    nrow(TD$event_data) # Only six audio events!

We will play with various detection thresholds: end_thr, SNR_thr, angle_thr, KPE and KME parameters. Try to reach 65 spectrograms extracted with a contour (i.e., Kalman curve) that best matches the audio event (answer below). The FFT size will be set at 512 samples.

    CATBIC <- read_audio("data/Catharus-bicknelli-54864.mp3")
    TD <- threshold_detection(
      CATBIC, threshold = 12, time_exp = 1, min_dur = 140, max_dur = 440, TBE = 10,
      EDG = 0.996, LPF = 10000, HPF = 1000, FFT_size = 512, FFT_overlap = 0.875,
      start_thr = 25, end_thr = 30, SNR_thr = 10, angle_thr = 125, duration_thr = 440,
      NWS = 1000, KPE = 1e-05, KME = 1e-05, settings = FALSE, acoustic_feat = TRUE,
      metadata = FALSE
    )

Let's take a look at the extracted audio features. Note that all the features are described and explained in the package vignette (`vignette("bioacoustics")`).

    # Acoustic features are stored in a data frame called event_data,
    # stored by order of detection.

    View(TD$event_data) # Contains the filename and the time of detection in the
                        # recording, and 26 extracted features.

The location (in number of samples) of the audio event in the recording is saved in a list.

    # Start and end of the 5th extracted audio event (in samples)
    c(TD$event_start[[5]], TD$event_end[[5]])

    # Remember you just have to divide by the sample rate to retrieve the time (s)
    c(TD$event_start[[5]], TD$event_end[[5]]) / slot(CATBIC, "samp.rate")

The amplitude (dB) and frequency (Hz) tracks (or bins) are also saved in a list. These can be used to build your own acoustic features.

    par(mar = c(1,1, 1, 1), oma = c(1, 1, 1, 1))

    # Amplitude track of the 5th audio event
    plot(TD$amp_track[[5]], type = "l")

    # Frequency track of the 5th audio event
    plot(TD$freq_track[[5]], type = "l")

The whole energy and frequency content can also be used to classify audio events instead of using acoustic features that may result in a loss of information. We will get there soon, but first, let's discover another detection function, here applied on echolocation calls of bats.

### Blob detection

The `blob_detection()` function will be used on a recording containing 10 bat echolocation calls from the Myotis genus. This function combines several image processing, filtering and image feature extraction. A blur and contrast boost is applied after mean background subtraction to increase the SNR of the audio event. The blob detection algorithm is applied on the processed spectrogram to detect the ROI (i.e., each preprocessed audio event). The blob detector simultaneously labels the connected FFT values and their contours in the spectrogram. Labelling is done in a single pass over the spectrogram, while contour points are revisited more than once and up to four times (see Chang et al., 2004). We will play with the arguments of this function to extract bat echolocation calls.

    # Access the arguments of the blob_detection function
    ?blob_detection

    # Use the bat recording stored in the package
    data(myotis)

    # Set each argument according to the targeted audio events
    BD <- blob_detection(
      myotis, # Either a path to an audio file (see ?read_audio), or a Wave object
      time_exp = 10, # Time expansion factor of 10 for time expanded recordings.
      min_dur = 1.5, # Minimum duration threshold of 1.5 milliseconds (ms)
      max_dur = 80, # Maximum duration threshold of 80 ms
      min_area = 40, # minimum number of 40 pixels in the blob
      TBE = 20, # Minimum time window between two audio events of 20 milliseconds
      EDG = 0.996, # Temporal masking with Exponential Decay Gain from 0 to 1
      LPF = slot(myotis, "samp.rate") * 10 / 2, # Low-Pass Filter at the Nyquist frequency
      HPF = 16000, # High-Pass Filter of 16 kHz
      FFT_size = 256, # Size of the Fast Fourrier Transform (FFT) window
      FFT_overlap = 0.875, # Percentage of overlap between two FFT windows

      blur = 2, # Gaussian smoothing function with a factor of 2 for blurring the spectrogram
      bg_substract = 20, # Foreground extraction with a mean filter applied on the spectrogram
      contrast_boost = 20, # Edge contrast enhancement filter of the spectrogram contour

      settings = FALSE, #  Save on a list the above parameters set with this function
      acoustic_feat = TRUE, # Extracts the acoustic and signal quality parameters
      metadata = FALSE, # Extracts on a list the metadata embedded with the Wave file
      spectro_dir = "./Spectros", # HTML page with spectrograms by order of detection
      time_scale = 0.1, # Time resolution of 2 ms for spectrogram display
      ticks = TRUE # Tick marks and their intervals are drawn on y-axis (frequencies)
    )

    # Get the number of extracted audio events
    nrow(BD$event_data)

Do not close the HTML page and tune the FFT size at 512. Let's play with the blur, contrast boost and background subtraction parameters to retrieve a number of 10 extracted echolocation calls.

    # Let's try various settings, starting with 512 FFT size instead of 256.
    BD <- blob_detection(
      myotis, time_exp = 10, FFT_size = 512, settings = FALSE, acoustic_feat = TRUE,
      metadata = FALSE, spectro_dir = "./Spectros", time_scale = 0.1, ticks = TRUE
    )

    # Take a look at the spectrograms and compare them with the previous extraction.
    nrow(BD$event_data) # Only 7 audio events!

Let's take a look at the extracted audio features. All the features are described and explained in the package vignette.

    # Acoustic features
    head(BD$event_data)

This data frame is, for now, the only available set of acoustic features with the `blob_detection()` function. However, it combines well with the `fspec()` to make image analysis.

Now that we have played with both detection functions with bird and bat vocalizations, let's go back to birds to explore batch analysis (i.e., with several recordings) and audio event classification.

## Batch analysis and classification

In this section, we will learn how to analyze several recordings at the same time and train a simple classifier (with training set) that will be used to classify new data (i.e., the test set).

We will work with 27 recordings of Catharus-bicknelli (n = 9), Passerella iliaca (n = 9), and Setophaga-magnolia (n = 9). We will split the extracted audio events in a 70 % training set (called “Train”) and 30 % test set (called “Test”).

Our target audio events are calls of Catharus-bicknelli. We will use the threshold detector previously configured for this species (see section 4.1.1).

    # Make sure you have the right path to the 27 MP3 recordings
    getwd()

    # Get the filepath for each MP3 file
    files <- dir("data", recursive = TRUE, full.names = TRUE, pattern = "[.]mp3$")

    # Detect and extract audio events
    TDs <- setNames(
      lapply(
        files,
        threshold_detection,
        threshold = 12, min_dur = 140, max_dur = 440, TBE = 50, LPF = 8000,
        HPF = 1500, FFT_size = 512, start_thr = 30, end_thr = 20, SNR_thr = 10,
        angle_thr = 125, duration_thr = 400, spectro_dir = NULL, NWS = 2000,
        KPE = 0.00001, time_scale = 2, EDG = 0.996
      ),
      basename(file_path_sans_ext(files))
    )

    # Keep only files with data in it
    TDs <- TDs[lapply(TDs, length) > 0]

    # Keep the extracted feature and merge in a single data frame for further analysis
    Event_data <- do.call("rbind", c(lapply(TDs, "[[", "event_data"), list(stringsAsFactors = FALSE)))
    nrow(Event_data) # 355 audio events extracted

    # Compute the number of extracted CATBIC calls
    sum(startsWith(Event_data$filename, "Cat"))

    # Add a "Class" column: "CATBIC" vs. other species of birds "OTHERS"
    classes <- as.factor(ifelse(startsWith(Event_data$filename, "Cat"), "CATBIC", "OTHERS"))
    Event_data <- cbind(data.frame(Class = classes), Event_data)

    # Get rid of the filename and time in the recording
    Event_data$filename <- Event_data$starting_time <- NULL

We now have the necessary dataset to train a classifier: we will train a Random Forest on the training set and validate the results on the test set.

    # Split the data in 60% Training / 40% Test sets
    train <- sample(1:nrow(Event_data), round(nrow(Event_data) * .6))
    Train <- Event_data[train,]

    test <- setdiff(1:nrow(Event_data), train)
    Test <- Event_data[test,]

    # Train a random forest classifier
    set.seed(666)
    rf <- randomForest(Class ~ duration + freq_max_amp + freq_max + freq_min +
                               bandwidth + freq_start + freq_center + freq_end +
                               freq_knee + fc + freq_bw_knee_fc + bin_max_amp +
                               pc_freq_max_amp + pc_freq_max + pc_freq_min +
                               pc_knee + temp_bw_knee_fc + slope + kalman_slope +
                               curve_neg + curve_pos_start + curve_pos_end +
                               mid_offset + smoothness + snr + hd + smoothness,
                       data = Train, importance = FALSE, proximity = FALSE,
                       replace = TRUE, ntree = 4000, mtry = 4)

    # Look at the confusion matrix of the training set
    rf$confusion # looks good, but...

    # Let's make predictions with our classifier on a test set
    table(Test[,1], predict(rf, Test[,-1], type = "response")) # not bad!

    # To look at the predictions
    head(predict(rf, Test[,-1], type = "prob"))

We are now able to use this simple, but proven robust, classifier to detect new calls of your target species.

# Deep learning classification with the R interface to Keras

We will use Keras in R which requires to install several packages in Python Guidelines to install Keras properly in R are available here

Let's now explore a ConvNet approach available on Keras. We will follow the approach of Hatami et al. (2017) to analyze time series as images with 2D ConvNets. The difference is that we will only perform max pooling at the last layer before activation and add batch normalization with dropouts at each layer.

    # Run if keras is installed on your machine

    # Build the training set
    Y_train <- to_categorical(as.integer(Train[,1]) - 1) # One hot encoding

    # X as matrix
    X_train <- as.matrix(Train[,-1])

    # Build the test set
    Y_test <- to_categorical(as.integer(Test[,1]) - 1)
    Y_test <- Y_test[,-1]
    X_test <- as.matrix(Test[,-1])

    # Build the sequential model
    mod0 <- keras_model_sequential()
    mod0 %>%
      # Input shape layer = c(samples, rows, cols, channels)
      layer_reshape(input_shape=ncol(X_train),target_shape=c(1,1,ncol(X_train))) %>%
      # First conv 2d layer with 128 neurons, kernel size of 8 x 8 and stride of 1 x 1
      layer_conv_2d(128, c(8,8), c(1,1), padding='same') %>%
      layer_batch_normalization() %>%
      layer_activation("relu") %>%
      layer_dropout(0.2) %>%
      # Second conv 2d layer with 256 neurons, kernel size of 5 x 5 and stride of 1 x 1
      layer_conv_2d(256, c(5,5), c(1,1), padding='same') %>%
      layer_batch_normalization() %>%
      layer_activation("relu") %>%
      layer_dropout(0.2) %>%
      # Third conv 2d layer with 128 neurons, kernel size of 3 x 3 and stride of 1 x 1
      layer_conv_2d(128, c(3,3), c(1,1), padding='same') %>%
      layer_batch_normalization() %>%
      layer_activation("relu") %>%
      layer_dropout(0.2) %>%
      # Average pooling layer
      layer_global_average_pooling_2d() %>%
      # Activation output layer with 2 classes
      layer_dense(units = ncol(Y_train),  activation='softmax')

    # Model compile
    mod0 %>% compile(loss = 'categorical_crossentropy',
                     optimizer = "adam",
                     metrics = "categorical_accuracy")

    # Add a callback to reduce the learning rate when reaching the plateau
    reduce_lr <- callback_reduce_lr_on_plateau(monitor = 'loss', factor = 0.5,
                                               patience = 50, min_lr = 0.0001)
    # Start learning
    mod0 %>% fit(X_train, Y_train, batch_size = 32, epochs = 50,
                 validation_data = list(X_test, Y_test),
                 verbose = 1, callbacks = reduce_lr)

    # Score on the test set
    score <- mod0 %>% evaluate(X_test, Y_test, batch_size = 32)
    score

Let's work a bit with the output to build a confusion matrix and use the predict function on the test set.

    # Look at predictions and build a confusion matrix
    Pred <- as.factor(predict_classes(mod0, X_test, batch_size = 32, verbose = 1))
    table(Y_test[,2], Pred)

    # To look at the prediction values
    Prob <- round(predict_proba(mod0, X_test, batch_size = 32, verbose = 1), 2)

We obtained a val_loss < 0.2 and val_categorical_accuracy > 0.94 which is acceptable, but not better than the simplest RF approach we used in section 3.2. Using only 26 acoustic features as model inputs instead of the whole spectrogram content (energy and frequency distribution, and harmonics) probably reduced the performances of the CNN model.

This tutorial is now complete. Comments and feedback are welcome:

Francois: francois.fabianek@wavx.ca
Jean: jean.marchal@wavx.ca
www.wavx.ca