---
title: "Detection, Extraction and Classification of Bird and Bat Vocalizations in R"
author: "Francois Fabianek, Jean Marchal"
date: "2018 October 30th"
output:
  rmarkdown::html_vignette:
    fig_caption: yes
    number_sections: yes
    toc: yes
    toc_depth: 4
  rmarkdown::pdf_document:
    fig_caption: yes
    number_sections: yes
    toc: yes
    toc_depth: 4
vignette: >
  %\VignetteIndexEntry{Tutorial: Detection, Extraction and Classification of Bird and Bat Vocalizations in R} 
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::knitr} 
---

_______
_______

# Load the necessary packages

Make sure you have the latest version of the packages:

```{r install_packages, message=FALSE, eval=FALSE}
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")
```

```{r load_packages, message=FALSE, eval=FALSE}
# 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:

```{r xeno1, message=FALSE, eval=FALSE}
df1 = query_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,]
```
  
```{r xeno2, message=FALSE, results='hold', eval=FALSE}
df2 = query_xc(qword ='Setophaga magnolia type:song cnt:"Canada"', download = FALSE)
df2 = df2[df2$Quality=="A",]
df2 = df2[1:9,]
```
  
```{r xeno3, message=FALSE, results='hold', eval=FALSE}
df3 = query_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)
```
  
```{r xeno4, eval=FALSE}
# Visualize your data frame
View(df)

# We will work in the R temp directory
wd <- tempdir()

# Create a data directory if it does not exist
data_dir <- file.path(wd, "data")

if(!dir.exists(data_dir))
  dir.create(data_dir)

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

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*:

```{r read_audio, eval=FALSE}
CATBIC <- read_audio(file.path(data_dir, "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.

___________________

# Extract GUANO metadata

[GUANO](https://guano-md.org/) stands for the "Grand Unified Acoustic Notation Ontology" and means to be a universal, extensible, open metadata format for bat (ultrasonic) and non-bat (audible, infrasonic) acoustic recordings [more here](https://github.com/riggsd/guano-r). GUANO format is now embeded directly in the WAV files generated from most Pettersson, Wildlife Acoustics, Titley Scientific acoustic recorders. It is possible to extract the metadata and GUANO embedded in the WAV file by using the `metadata()` function.

```{r metadata,eval=FALSE}

metadata(CATBIC)

```
<br>


Now that we have explored a WAV file, we will use the Fast Fourrier Transform (FFT) to compute a frequency-time representation of the recording, called a spectrogram. A spectrogram being the representation of the spectrum of frequencies in a recording as they vary through time. This representation, although not optimal, is still commonly 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.

```{r spectro0, eval=FALSE}
# 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:  
  
```{r help, eval=FALSE}
# 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.

```{r spectro1, eval=FALSE}
# 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:

```{r filter, eval=FALSE}
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*).  

```{r threshold_help, eval=FALSE}
# Access the arguments of the threshold_detection function
?threshold_detection
```

```{r threshold1, eval=FALSE}
# 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
  min_TBE = 10, # Minimum time window between two audio events of 10 milliseconds
  max_TBE = 5000, # Maximum time window between two audio events, here 5 seconds
  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 = file.path(tempdir(), "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$data$event_data)
```

Let the HTML page open with the 57 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.

```{r threshold2, eval=FALSE}
# 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, 
  min_TBE = 10, max_TBE = 5000, 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 = file.path(tempdir(), "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$data$event_data) # Only three audio events!
```

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

```{r threshold3, eval=FALSE}
CATBIC <- read_audio(file.path(data_dir, "Catharus-bicknelli-54864.mp3"))
TD <- threshold_detection(
  CATBIC, threshold = 12, time_exp = 1, min_dur = 140, max_dur = 440, min_TBE = 10, 
  max_TBE = Inf, EDG = 0.996, LPF = 10000, HPF = 1000, FFT_size = 256, FFT_overlap = 0.875, 
  start_thr = 22, 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")`).

```{r features1, eval=FALSE}
# Acoustic features are stored in a data frame called event_data,
# stored by order of detection.

View(TD$data$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. 

```{r features2, eval=FALSE}
# Start and end of the 5th extracted audio event (in samples)
c(TD$data$event_start[[5]], TD$data$event_end[[5]])

# Remember you just have to divide by the sample rate to retrieve the time (s)
c(TD$data$event_start[[5]], TD$data$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.

```{r features3, eval=FALSE}
par(mar = c(1,1, 1, 1), oma = c(1, 1, 1, 1))

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

# Frequency track of the 5th audio event
plot(TD$data$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](https://doi.org/10.1016/j.cviu.2003.09.002). We will play with the arguments of this function to extract bat echolocation calls.

```{r blob0, eval=FALSE}
# Access the arguments of the blob_detection function
?blob_detection
```

```{r blob1, eval=FALSE}
# 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
  min_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 = file.path(tempdir(), "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$data$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.

```{r blob2, eval=FALSE}
# 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 = file.path(tempdir(), "Spectros"), time_scale = 0.1, ticks = TRUE
) 

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

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

```{r blobfeat1, eval=FALSE}
# Acoustic features
head(BD$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).

```{r classification1, eval=FALSE}
# Get the filepath for each MP3 file
files <- dir(data_dir, 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, min_TBE = 50, max_TBE = Inf,
    LPF = 8000, HPF = 1500, FFT_size = 256, 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, function(x) length(x$data)) > 0]

# Keep the extracted feature and merge in a single data frame for further analysis
Event_data <- do.call("rbind", c(lapply(TDs, function(x) x$data$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.

```{r classification2, eval=FALSE}
# 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](https://www.python.org/downloads/)
Guidelines to install Keras properly in R are available [here](https://keras.rstudio.com/)

Let's now explore a ConvNet approach available on Keras. We will follow the approach of [Hatami et al. (2017)](https://arxiv.org/pdf/1710.00886.pdf) 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.

```{r keras1, eval=FALSE}
# Run if keras is installed on your machine
library(keras)

# 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.

```{r keras2, eval=FALSE}
# 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](https://www.wavx.ca)


_______
_______