inst/doc/sonar.R

In bioSNR: Bioacoustic Basic Operations with Decibels and the Passive Sonar Equation

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup, eval=FALSE--------------------------------------------------------
#  #Stable - Install package from CRAN
#  install.packages("bioSNR")
#  
#  #Unstable - Install package from Github repository
#  devtools::install_github("MattyD797/bioSNR")
#  
#  #Attach package namespace to active libraries in Rstudio
#  library(bioSNR)

## ----setupR, eval=TRUE, include=FALSE-----------------------------------------
library(bioSNR)

## ----wenz, echo=FALSE, fig.align = 'center', message=FALSE--------------------
freqBand <- c(0,0)
shipT <- -1
seaState <- -1
wSpeed <- 0
boolR <- T
    # This code re-traces the empirical Wenz curves, that describe Oceanic ambient noise as a function of frequency.
    freqMax <- 10000
    freq <- pracma::logspace(0,6,freqMax)


    #Geophysical contribution for Freq<10 Hz
    #N_geophys = 107 - 30log10(f)
    intervalG <- freq[1:tail(which(freq<10)+1, 1)]
    NL_geophys <- 107 - 30 * log10(intervalG)

    NL_geophys2 <- cbind(intervalG, NL_geophys)
    colnames(NL_geophys2)[1] <- "freq"

    ## Contribution of ship traffic (5 <= f <= 500 Hz)
    # - 1 - 2 low ship traffic
    # - 3-4-5 standard ship traffic
    # - 6 - 7 heavy ship traffic
    # - 8 - 9 intense ship traffic

    # N_ship = 75 -40log10(f/30) + 5(n_shiptraffic - 4)

    shipLines <- function(trafficIndex, freqI=freq, freqmax = 10000){

      interval <- freqI[head(which(freqI > 5), 1):tail(which(freqI < 500), 1)]
      NL_ship<- 76 - (20*(log10(interval/30))^2) + 5*(trafficIndex-4)
      return(NL_ship)
    }

    shipTraffic <- sapply(c(1:9), shipLines)
    shipTraffic2 <- cbind(freq[head(which(freq > 5), 1):tail(which(freq < 500), 1)], shipTraffic)
    colnames(shipTraffic2) <- c("freq",
                                "shipTraffic.1",
                                "shipTraffic.2",
                                "shipTraffic.3",
                                "shipTraffic.4",
                                "shipTraffic.5",
                                "shipTraffic.6",
                                "shipTraffic.7",
                                "shipTraffic.8",
                                "shipTraffic.9")

    ## Wind speed and sea state  (80 <= f <= 300 kHz)
    # Wind speed in knots - related to sea state as Vkn = 5 * sea state

    windSpeed <- function(windSpeed, freqI=freq, freqmax = 10000){
      #between 80 - 1000 Hz
      # < 1000 Hz
      interval1 <- head(which(freqI > 80)-1, 1):(tail(which(freqI < 1000)-1, 1))
      NL_wind <- 44 + sqrt(21*windSpeed) + 17 * (3 - log10(freqI[interval1])) * (log10(freqI[interval1])-2)
      # >= 1000 Hz
      interval2 <- head(which(freqI > 1000)-1, 1):tail(which(freqI < 300000)-1, 1)
      NL_wind <- append(NL_wind, 95 + sqrt(21*windSpeed) - 17*log10(freqI[interval2]))
      return(NL_wind)
    }

    if(seaState != -1){
      wSpeed = seaState * 5
    }

    wind <- sapply(c(0, 10, 15, 20, 25, 30, 40, 45), windSpeed)
    intervalW <- freq[head(which(freq > 80)-1, 1):tail(which(freq < 300000)-1, 1)]
    wind2 <- cbind(intervalW, wind)
    colnames(wind2) <- c("freq", "seaState.0", "seaState.10", "seaState.15", "seaState.20", "seaState.25", "seaState.30", "seaState.40", "seaState.45")

    ## Thermal noise (f > 50 kHz)
    # NLthermal = -75 + 20log10(f)
    intervalT <- freq[head(which(freq > 50000), 1):freqMax]
    NL_thermal <- -75 + 20 * log10(intervalT)
    NL_thermal2 <- cbind(intervalT, NL_thermal)
    colnames(NL_thermal2)[1] <- "freq"

    freq <- as.data.frame(freq)
    colnames(freq) <- "freq"

    #Solve for given bandwidth
    bandW <- freq[freq>=freqBand[1] & freq <= freqBand[2]]
    gpInputF <- vector()
    stInputF <- vector()
    wInputF <- vector()
    tnInputF <- vector()
    gpInput <- vector()
    stInput <- vector()
    wInput <- vector()
    tnInput <- vector()

    for(f in bandW){
      #Only geophysical SL
      if(f < 10){
        gpInput<- c(gpInput, 107 - 30 * log10(f))
        gpInputF <- c(gpInputF, f)

      }
      if (f >= 5 & f <= 500){
        if(shipT < 0 ){
          if(f <= 10){
            stop("ERROR: shipT is not specified")
          }
        } else {
          stInput <- c(stInput, 76 - (20*(log10(f/30))^2) + 5*(shipT-4))
          stInputF <- c(stInputF, f)
        }
      }
      if(f>=80 & f < 1000){
        wInput <- c(wInput, 44 + sqrt(21*wSpeed) + 17 * (3 - log10(f)) * (log10(f)-2))
        wInputF <- c(wInputF, f)
      }
      if(f>=1000 & f <= 300000){
        wInput <- c(wInput, 95 + sqrt(21*wSpeed) - 17*log10(f))
        wInputF <- c(wInputF, f)
      }
      if(f>50000){
        tnInput <- c(tnInput,-75 + 20 * log10(f))
        tnInputF <- c(tnInputF, f)
      }
    }

    gpInput2 <- cbind(gpInputF, gpInput)
    colnames(gpInput2)[1] <- "freq"
    stInput2 <- cbind(stInputF, stInput)
    colnames(stInput2)[1] <- "freq"
    wInput2 <- cbind(wInputF, wInput)
    colnames(wInput2)[1] <- "freq"
    tnInput2 <- cbind(tnInputF, tnInput)
    colnames(tnInput2)[1] <- "freq"


    df_list <- list(freq, NL_geophys2, shipTraffic2, wind2, NL_thermal2, gpInput2, stInput2, wInput2, tnInput2)
    dt <- Reduce(function(x, y) merge(x, y, by = "freq", all = TRUE), df_list)

    dt <- dplyr::mutate(dt, gpInput = as.numeric(gpInput),
             stInput = as.numeric(stInput),
             wInput = as.numeric(wInput),
             tnInput = as.numeric(tnInput))

    suppressWarnings(print(ggplot2::ggplot(dt, ggplot2::aes(freq)) +
      ggplot2::geom_line(ggplot2::aes(y=shipTraffic.1))+
      ggplot2::geom_line(ggplot2::aes(y=shipTraffic.2))+
      ggplot2::geom_line(ggplot2::aes(y=shipTraffic.3))+
      ggplot2::geom_line(ggplot2::aes(y=shipTraffic.4))+
      ggplot2::geom_line(ggplot2::aes(y=shipTraffic.5))+
      ggplot2::geom_line(ggplot2::aes(y=shipTraffic.6))+
      ggplot2::geom_line(ggplot2::aes(y=shipTraffic.7))+
      ggplot2::geom_line(ggplot2::aes(y=shipTraffic.8))+
      ggplot2::geom_line(ggplot2::aes(y=shipTraffic.9))+
      ggplot2::scale_x_continuous(trans='log10', ,
                         breaks = scales::trans_breaks("log10", function(x) 10^x),
                         labels = scales::trans_format("log10", scales::math_format(10^.x)),
                         n.breaks = 25)+
      ggplot2::geom_line(ggplot2::aes(y=NL_geophys))+
      ggplot2::geom_line(ggplot2::aes(y=seaState.0))+
      ggplot2::geom_line(ggplot2::aes(y=seaState.10))+
      ggplot2::geom_line(ggplot2::aes(y=seaState.15))+
      ggplot2::geom_line(ggplot2::aes(y=seaState.20))+
      ggplot2::geom_line(ggplot2::aes(y=seaState.25))+
      ggplot2::geom_line(ggplot2::aes(y=seaState.30))+
      ggplot2::geom_line(ggplot2::aes(y=seaState.40))+
      ggplot2::geom_line(ggplot2::aes(y=seaState.45))+
      ggplot2::geom_line(ggplot2::aes(y=NL_thermal))+
      ggplot2::ylab(as.expression(bquote("Spectrum Level (dB ref 1 \u00B5" ~ Pa^2 ~ ')')))+
      ggplot2::xlab("Frequency (Hz)") +
      ggplot2::geom_point(ggplot2::aes(y=gpInput))+
      ggplot2::geom_point(ggplot2::aes(y=stInput))+
      ggplot2::geom_point(ggplot2::aes(y=wInput))+
      ggplot2::geom_point(ggplot2::aes(y=tnInput))+
      ggplot2::geom_vline(xintercept = freqBand[1])+
      ggplot2::geom_vline(xintercept = freqBand[2])+
      ggplot2::annotation_logticks(sides = "b")))

## ----ex1, fig.align = 'center', warning=FALSE---------------------------------
specLvlGraph(c(28,33), ship=4,seaState = 1, wSpeed = 10, boolR = T)[[2]]

## ----ex2, fig.align = 'center'------------------------------------------------
#Source Level
SL <- 195
#Noise Level - from above
NL <- 82.9897
#Detection Threshold
DT <- 10
#Transition Range <- depth/2
TR <- 2500

rmax(sl=SL,nl=NL,dt=DT,d=5000, xaxis=10000000)