classifier.Rmd

---
title: "Classifying drone RF signals with statistical learning and a small data set"
author: "Ben Horvath"
date: "July 2022"
output:
  html_document:
    toc: true
    number_sections: true
    theme: null
    css: ./static/architect.css
    template: ./static/architect.html
    pandoc_args: [
      "--mathjax", "",
      "--variable", "mathjax-url:https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"]
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(warning=FALSE, message=FALSE)
options(scipen=999, digits=4)

setwd('~/dronerf_classifier/')

library(caret)
library(dplyr)
library(ggplot2)
library(knitr)
library(stringr)
library(tuneR)
library(seewave)
library(png)

DATA_PATH <- './data'

dir.create(DATA_PATH, showWarnings=FALSE)

IMG_SIZE <- 122
```

This brief note puts together a couple (non-deep learning) algorithms to classify RF signals using a small open-source data set.This work agrees with Medaiyese, et al. (2021) that large labeled data sets and complicated deep learning may not be essential for classifying drone RF signals.

GitHub repo: [http://github.com/benhorvath/dronerf_classifier](http://github.com/benhorvath/dronerf_classifier/)



# Data and summary

The data is available for download at: [`DroneRF` dataset: A dataset of drones for RF-based detection, classification, and identification](https://data.mendeley.com/datasets/f4c2b4n755/1). Uncompressed, it is 43 gigabytes in size. It is summarized in Allahham, et al. (2019):

> ... the `DroneRF` dataset: a radio frequency (RF) based dataset of drones functioning in different modes, including off, on and connected, hovering, flying, and video recording. The dataset contains recordings of RF activities, composed of 227 recorded segments collected from 3 different drones, as well as recordings of background RF activities with no drones. The data has been collected by RF receivers that intercepts the drone's communications with the flight control module. The receivers are connected to two laptops, via PCIe cables, that runs a program responsible for fetching, processing and storing the sensed RF data in a database. 

The data set is very clean, collected under laboratory conditions using expensive scientific equipment. The entire 2.4 Ghz wifi bandwidth was monitored with two RF receivers, producing both a lower and higher half of the frequency band. Medaiyese, et al. (2021) determined the lower alone contained sufficient information for good classification performance. Each of the 227 records is equal to or less than 10 seconds of recording, indicating this data set is fairly small.

After downloading and decompressing the data, the first step is to process the records and transform them into spectrograms of size $122 x 122$ pixels. Each spectrograms' pixel values (intensity) is taken, and then converted into a $1 x (122*122) = 1 x 14,884$ matrix. This produces a feature  matrix appropriate for input to algorithms. The data is split into a training partition and a final test partition. R's `caret` library is used to train several algorithms amenable to $n << p$ classification tasks with $5-$fold cross-validation. These are GLM with ElasticNet regularization (Zou and Hastie 2005; Friedman, Hastie, and Tibshirani 2010) and Random Forest. Performance is quite satisfactory on both binary and multi-class tasks, especially so considered the limited hardware thrown at the problem and the size of the data set. This work agrees with Medaiyese, et al. (2021) that large labeled data sets and complicated deep learning may not be essential for classifying drone RF signals.

Three models are developed and tested on a reserved hold-out set:

* $M_0$: Binary task, binomial GLM with ElasticNet regularization (R library: `glmnet`); hold-out performance:
  - precision: 0.67
  - recall: 1.0
  - F-score: 0.80
  - balanced accuracy: 0.95
* $M_1$: Binary task, Random Forest (R library: `ranger`); hold-out performance:
  - precision: 0.91
  - recall: 1.0
  - F-score: 0.95
  - balanced accuracy: 0.99
* $M_2$: Multiclass task, Random Forest; hold-out performance:
  - mean balanced accuracy: 0.88
  - mean sensitivity: 0.83
  - mean specificity: 0.92

Model performance could easily be improved with slightly more powerful hardware, which would allow spectograms to encode more information, i.e., larger than $122 x 122$ pixels.




# Load and process data

The process here is to 'summarize' each ~100 megabyte RF recording as a ~4 kilobyte spectrogram PNG image, and then the spectrogram into a ~100 kilobyte $1 x (122*122)$ matrix.

Create necessary directories (here input refers to input to classifier):

```{r}
dronerf_path <- file.path(DATA_PATH, 'dronerf')
dir.create(dronerf_path, showWarnings=FALSE)

input_path <- file.path(DATA_PATH, 'input')
dir.create(input_path, showWarnings=FALSE)

input_files <- file.path(dronerf_path, list.files(dronerf_path, recursive=TRUE))

lower <- input_files[str_detect(input_files, '_L[\\d]?/') == TRUE]

sample_rate <- 40 * 10^6
```



## Create spectrograms

Process lower half and save as PNG to disk:

```{r, eval=FALSE}
name_map <- list('AR drone' = 'ar',
                 'Bepop drone' = 'bepop',
                 'Background RF activites' = 'bg',
                 'Phantom drone' = 'phantom')
                 
for (f in lower) {
  
  f_prefix <- str_extract(f, 'AR drone|Bepop drone|Background RF activites|Phantom drone')  # NOTE the typo in 'activites'
  drone <- name_map[[ f_prefix ]]
  
  f_suffix <- str_match_all(f, '(\\d{5}L_\\d+).csv$')[[1]][2]
  
  sig <- scan(f, what=character(1), sep=',') %>%
    as.numeric %>%
    tuneR::Wave(samp.rate=sample_rate) %>%
    tuneR::normalize(unit = c('1'))
  
  y = sig@left
  
  N <- 32e-3*sample_rate
  
  # ccompute spectorgram
  p <- seewave::ggspectro(y, sample_rate, wl=N, wn='hamming', ovlp=50, fftw=TRUE) +
    geom_tile(aes(fill = amplitude))+
    scale_fill_gradient(low='white', high='black') +
    theme(axis.line=element_blank(),
          axis.text.x=element_blank(),
          axis.text.y=element_blank(),
          axis.ticks=element_blank(),
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          legend.position='none',
          panel.background=element_blank(),
          panel.border=element_blank(),
          panel.grid.major=element_blank(),
          panel.grid.minor=element_blank(),
          plot.background=element_blank())
  
  export_path <- file.path(input_path, drone)
  dir.create(export_path, showWarnings=FALSE)
  
  export_file <- paste(file.path(export_path, f_suffix), '.png', sep='')
  
  ggsave(export_file, p, width=IMG_SIZE, height=IMG_SIZE, units='px')
  
}
```



## Spectrograms to feature matrix

Convert each PNG to a matrix, then combine to form a feature matrix $X$ for machine learning:

```{r}
standardize <- function(x) (x - min(x)) / (max(x) - min(x))

X_list <- list()
y <- list()

for (drone in c('ar', 'bepop', 'bg', 'phantom')) {
  
  pngs_path <- file.path(input_path, drone)
  drone_pngs_path <- file.path(pngs_path, list.files(pngs_path))
  
  for (png_path in drone_pngs_path) {
    
    img <- readPNG(png_path)
    x <- img[,,1]
    x <- t(apply(x, 2, rev))  # rotate: not necessary except for plotting correct orientation
    x <- standardize(x)
    dim(x) <- c(1, IMG_SIZE*IMG_SIZE)
    
    X_list[[1+length(X_list)]] <- x
    y[[1+length(y)]] <- drone
    
  }
  
}

rm(x, sig, p); gc()  # free up some memory

# Combine into single matrix for input to algorithm
X <- do.call(rbind, X_list) %>% data.frame
colnames(X) <- paste('p', 1:(IMG_SIZE*IMG_SIZE), sep='')

rm(X_list); gc()
```



# Classifier: Drone or no drone

## Build test and training sets

We will use $5-$fold cross-validation to train and tune the learning algorithms, but retain 25% of the data set for final testing.

For all models, under-sampling will be used to balance the classes. 

```{r}
df <- X %>%
  mutate(y=unlist(y)) %>%
  select(y, everything())

df <- df %>%
  mutate(y = as.factor(if_else(y == 'bg', 'bg', 'drone')))

set.seed(1804)
train_ix <- createDataPartition(df$y, p=0.75, list=FALSE) %>% as.numeric

train <- df[train_ix,]
test <- df[-train_ix,]

# check distributions are similar
prop.table(table(df$y))
prop.table(table(df$y[train_ix]))
prop.table(table(df$y[-train_ix]))
```



## $M_0$: Binomial GLM with Elasticnet regularization

Preprocessing _occurs within each fold_, and includes: Eliminating near-zero-variance predictors, scaling and centering:

```{r}
Sys.time()

cv5 <- trainControl(method='repeatedcv',
                     number=5,
                     classProbs=TRUE,
                     sampling='down',
                     summaryFunction = twoClassSummary,
                     savePredictions=TRUE)

m0_grid <- expand.grid(alpha = c(1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 0.1, 1),
                       lambda = c((1:5) / 10))

m0 <- train(y ~ .,
            preProc = c('nzv', 'center', 'scale'),
            data=train, 
            tuneGrid=m0_grid,
            method='glmnet',
            trControl=cv5,
            metric='ROC')

Sys.time()

m0
```

 1e-04  0.5     0.9928571  1.0000000  0.9142857

## $M_1$: Random forest

```{r}
Sys.time()

m1_cv <- trainControl(method='repeatedcv',
                     number=5,
                     classProbs=TRUE,
                     sampling='down',
                     summaryFunction = twoClassSummary,
                     savePredictions=TRUE)

# common mtry values are sqrt(p) and log2(p)
m1_tune <- expand.grid(mtry=c(14, 122), splitrule=c('gini'), min.node.size=c(5, 10))

m1 <- train(y ~ .,
            data=train,
            tuneGrid=m1_tune,
            method='ranger',
            num.trees=1000,
            trControl=m1_cv,
            metric='ROC')

Sys.time()

m1
```

## Evaluating $M_0$ and $M_1$ on hold out set

Performance metrics for both algorithms on the same hold-out set is presented below. Both perform similarly, with a slight edge to the Random Forest model, which also trains moderately faster. 


```{r}
m0_pred <- predict(m0, test)
m1_pred <- predict(m1, test)

confusionMatrix(m0_pred, test$y)
confusionMatrix(m1_pred, test$y)


```


# Classifier: Drone model + background


## Build test and training sets

```{r}
df <- X %>%
  mutate(y=unlist(y)) %>%
  select(y, everything())

df <- df %>%
  mutate(y = as.factor(y))

set.seed(1804)
train_ix <- createDataPartition(df$y, p=0.75, list=FALSE) %>% as.numeric

train <- df[train_ix,]
test <- df[-train_ix,]
```


## $M_2$: Random forest: All drones v. background

```{r}
Sys.time()

m2_cv <- trainControl(method='repeatedcv',
                     number=5,
                     classProbs=TRUE,
                     sampling='down',
                     summaryFunction = multiClassSummary,
                     savePredictions=TRUE)

# common mtry values are sqrt(p) and log2(p)
m2_tune <- expand.grid(mtry=c(14, 122), splitrule=c('gini'), min.node.size=c(5, 10))

m2 <- train(y ~ .,
            data=train,
            tuneGrid=m2_tune,
            method='ranger',
            num.trees=1000,
            trControl=m2_cv,
            metric='AUC')

Sys.time()

m2
```

## Evaluating $M_2$ on hold out set

$M_2$ performs less successfully than the previous models on binary target, but still does reasonably well. 

```{r}
m2_pred <- predict(m2, test)

confusionMatrix(m2_pred, test$y)
```



# References

* Allahham, MHD Saria, Mohammad F. Al-Sa'd, Abdulla Al-Ali, Amr Mohamed, Tamer Khattab, and Aiman Erbad. 2019. '`DroneRF` dataset: A dataset of drones for RF-based detection, classification and identification.' _Data in Brief_ 26: 104313.

* Friedman, Jerome, Trevor Hastie, and Rob Tibshirani. 2010. 'Regularization paths for generalized linear models via coordinate descent.' _Journal of Statistical Software_ 33, no. 1: 1.

* Medaiyese, Olusiji O., Abbas Syed, and Adrian P. Lauf. 2021. 'Machine learning framework for RF-based drone detection and identification system.' _2021 2nd International Conference on Smart Cities, Automation \& Intelligent Computing Systems_. IEEE.

* Zou, Hui, and Trevor Hastie. 2005. 'Regularization and variable selection via the elastic net.' _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 67, no. 2: 301--20.