###########################################################################################
##    ANALYSIS PER CALL
#####################################################################################

## Notes:
# same among all versions 0-4, except for a set of square brackets at "filPeaNot =", probably not important
# NB: check you use version 1.8.2.tar.gz	2020-10-17 17:30	3.9M for
# 'soundgen' package! it changes a lot among versions and behaves inconsistently

#############################################################
##Packages:
require(soundgen)
require(seewave)
require(tuneR)

###################################################################################
## CHANGE FILE-SPECIFIC PARAMETERS IN FOLLOWING SECTION:

## Set the location of the files to be analysed and make list of the files, below:
setwd("")

## new list of file names, only .wav / .WAV
fileName = list.files()
fileName = fileName[(grepl(".WAV",fileName))|(grepl(".wav",fileName))]

## recording ID: file name without _# or extension
recordID = lapply(fileName, function(x){
  gsub("_.*", "", x)
})

## call Nr: # after the '_' in the file name (probably not the most efficient way)
calNr = lapply(fileName, function(x){
  gsub(".*_", "", x)
})
# remove extension
calNr = lapply(calNr, function(x){
  gsub("(.*)\\..*$", "\\1", x)
})
# make numeric
calNr = as.numeric(unlist(calNr))

## call ID: file name without extension
calID = lapply(fileName, function(x){
  gsub("(.*)\\..*$", "\\1", x)
})

###########################################################################################
##  convert Wav files to Wave objects (TuneR)

waves = lapply(fileName, function(x){
  readWave(x)
})

###########################################################################################
## Frequency analysis of whole call
###########################################################################################

## list of peaks for power spectrum for each l
frePea = lapply(waves, function(x){
  fpeaks(meanspec(x, wl = 1024, ovlp = 50, plot =  F),
         plot = F)
})

## dominant freqcy 
domFre = lapply(frePea, function(x){
  x[,1][x[,2] == max(x[,2])]
})

## dom freq amplitude
domFreAmp = lapply(frePea, function(x){
  x[,2][x[,2] == max(x[,2])]
})

## basal frequency: won't be below 300 Hz
basFre = lapply(frePea, function(x){
  min(x[,1][x[,1] >= 0.3])
})

# basal freq amp
basFreAmp = lapply(frePea, function(x){
  min(x[,2][x[,1] >= 0.3])
})

# basal freq amplitude % reduction from domFreq
basFreAmpRed = lapply(1:length(fileName), function(x){
  (domFreAmp[[x]]-basFreAmp[[x]]) / domFreAmp[[x]]*100
})

## lower 90% bandwidth peak
banWidLow = lapply(frePea, function(x){
  min(x[,1][
    x[,2] >= (0.1*(x[,2][x[,2] == max(x[,2])]))
  ]) 
})

## upper 90% bandwidth peak
banWidUpp = lapply(frePea, function(x){
  max(x[,1][
    x[,2] >= (0.1*(x[,2][x[,2] == max(x[,2])]))
  ]) 
})


#################################################################
## 2nd emphasised freqcy: 
## if > 1kHz different from dom freq and < 6 kHz
## (and its amplitude):

# filter step
filPea = lapply(1:length(fileName), function(x){
  frePea[[x]][
    frePea[[x]][,1] < 6 &
      frePea[[x]][,1] > (domFre[[x]]+1) | 
      frePea[[x]][,1] < (domFre[[x]]-1)
    ,,drop = F]
})

# 2nd emphasised freq
empFre2 = lapply(filPea, function(x){
  x[,1][x[,2] == max(x[,2])]
})

# 2nd emph freq amplitude
empFre2Amp = lapply(filPea, function(x){
  x[,2][x[,2] == max(x[,2])]
})

# 2nd freq amplitude % reduction from domFreq
empFre2AmpRed = lapply(1:length(fileName), function(x){
  (domFreAmp[[x]]-empFre2Amp[[x]]) / domFreAmp[[x]]*100
})

###########################################################################################
## Segment call into notes for further analysis
###########################################################################################


###########################################################################################
## Filer out noise and ready the calls for segmentation

## filter out frequencies below 0.9 kHz (mostly noise) 
## converting wav to vector
wavesF = lapply(waves, function(x){
  ffilter(x, from = 0, to = 900, bandpass = F, 
          wl = 1024, ovlp = 50,  
          output = "matrix")
})

##  calculate the amplitudes for the frequency-filtered noise for time
amplts = lapply(1:length(wavesF), function(x){
  env(wavesF[[x]], f = waves[[x]]@samp.rate, envt = "abs", msmooth = c(32, 50), 
      norm = F, plot = F)
})

###########################################################################################
##  set the syllable threshold value at a standard, normalised to the loudest 2%
##    of the call (rather than normalising to the mean amplitude of the whole sound,
##    which is the default for function "segment" below, and causes too much vairability)

# scaling factor: mean amplitude of top 1%
SF = lapply(amplts, function(x)
{mean(na.omit(tail(sort(x), (nrow(x)/100))))})


# syllable threshold for segmentation ("segment" below), normalising 
# mean amplitude to SF:
ST = lapply(1:(length(fileName)), function(x){
  SF[[x]]/20/mean(amplts[[x]], na.rm = T)
}) 


###########################################################################################
##  Segment the sound into syllables using normalised threshold value (ST), for 
##  initial note break-up times

seg = lapply(1:(length(fileName)), function(x){
  segment(wavesF[[x]][,1], samplingRate = waves[[x]]@samp.rate, 
          windowLength = (3=2/(waves[[x]]@samp.rate/1000)), 
          overlap = 50, shortestSyl = 1, 
          sylThres = ST[[x]], shortestPause = 15,# interburst = 50, 
          burstThres = 1, 
          plot = T, main = fileName[[x]]
  )})

## Filter the segments to only syllables (remove bursts) > 15 MS, 
## to cut out remaining noise bursts

seg = lapply(seg, function(x){
  x$syllables[(x$syllables$sylLen > 15),]
})


#################################################################################################################
## Frequency analyses per note
###########################################################################################


## At sampling rate of x kHz, a single sample lasts MS milliseconds, 
## hence the first sample is from 0 to MS milliseconds, the second from
## MS to MS*2 milliseconds, etc:
MS = lapply(1:length(fileName), function(x){
  1/(waves[[x]]@samp.rate/1000)
})


## list of peaks for spectograms for each note, within each call
frePeaNot = lapply(1:length(fileName), function(x){
  lapply(1:nrow(seg[[x]]), function(y){
    fpeaks(meanspec(
      waves[[x]][
        (if((seg[[x]]$start[y]/MS[[x]]-15/MS[[x]])<=0){
          1}else{(seg[[x]]$start[y]/MS[[x]]-15/MS[[x]])}): #the sample number, e.g., 4th sample, that starts at that particular time in milliseconds, (should actually be seg[[x]]$start[1])/MS+1), but we'll start ~15 milliseconds early and end 10 ms late just to give time for an adequate time analysis; built-in safety against starting before the start of the recording
          (seg[[x]]$end[y]/MS[[x]]+15/MS[[x]])],
      wl = 1024, ovlp = 50, plot =  F),
      plot = F)
  })
})

## dominant freqcy (per note per call)
domFreNot = lapply(1:length(fileName), function(x){
  sapply(1:nrow(seg[[x]]), function(y){
    frePeaNot[[x]][[y]][,1][
      frePeaNot[[x]][[y]][,2] == 
        max(frePeaNot[[x]][[y]][,2])
    ]
  })
})

## basal frequency: won't be below 350 Hz
basFreNot = lapply(1:length(fileName), function(x){
  lapply(1:nrow(seg[[x]]), function(y){
    min(frePeaNot[[x]][[y]][,1][
      frePeaNot[[x]][[y]][,1] >= 0.35
    ])
  })
})

#################################################################
## 2nd emphasised freqcy (per note per call): 
## if > 1kHz different from dom freq and < 6 kHz
## (and its amplitude):

# filter step
filPeaNot = lapply(1:length(fileName), function(x){
  lapply(1:nrow(seg[[x]]), function(y){
    frePeaNot[[x]][[y]][
      frePeaNot[[x]][[y]][,1] < 6 &
        frePeaNot[[x]][[y]][,1] > (domFreNot[[x]][[y]]+1) | 
        frePeaNot[[x]][[y]][,1] < (domFreNot[[x]][[y]]-1)
      ,,drop = F]
    
  })
})


# 2nd emphasised freq
empFre2Not = lapply(1:length(fileName), function(x){
  lapply(1:nrow(seg[[x]]), function(y){
    filPeaNot[[x]][[y]][,1][
      filPeaNot[[x]][[y]][,2] == max(filPeaNot[[x]][[y]][,2])
    ]
  })
})


##################################################################################################
## Condense per note frequency measures to reported variables
#############################################################################################

## median dominant frequency for all notes
domFreNotMd = lapply(domFreNot, function(x){
  median(x)})

## range of dominant frequency for all notes
domFreNotRan = lapply(domFreNot, function(x){
  max(x)-min(x)})

## first note dominant frequency deviance
domFre1Dev = lapply(domFreNot, function(x){
  x[1]-median(x)
})


###########################################################################################
##  further variable columns for analysis (some already done above)
###########################################################################################


###########################################################################################
##note/call timing variables:

#number of notes:
nNot = lapply(seg, function(x){nrow(x)})

#call duration
calDur = lapply(seg, function(x){
  (x$end[nrow(x)] - x$start[1])/1000
})

#note rate (/second):
notRat = lapply(1:(length(fileName)), function(x){
  (nNot[[x]]-1)/((seg[[x]]$start[nNot[[x]]]-seg[[x]]$start[1])/1000)})

#median note duration (milliseconds):
notDurMd = lapply(seg, function(x){median(x$sylLen)
})

#note duration range %
notDurRan = lapply(seg, function(x){
  (max(x$sylLen) - min(x$sylLen)) / median(x$sylLen) * 100
})

##first note duration deviance as percentage:
notDur1Dev = lapply(1:length(seg), function(x){
  (seg[[x]]$sylLen[1] - notDurMd[[x]]) / notDurMd[[x]]*100
})

#longest note number
lonNotNr = lapply(seg, function(x){
  mean(x$syllable[x$sylLen==max(x$sylLen)])})

#shortest note number
shoNotNr = lapply(seg, function(x){
  mean(x$syllable[x$sylLen==min(x$sylLen)])})


###########################################################################################
##inter-note timing analysis:

#median Inter-note interval (milliseconds):
intNotIntMd = lapply(seg, function(x){median(x$pauseLen, na.rm = T)})

#inter-note interval range percentage:
intNotIntRan = lapply(1:length(seg), function(x){
  (max(seg[[x]]$pauseLen, na.rm = T) - min(seg[[x]]$pauseLen, na.rm = T)) /
    (intNotIntMd[[x]])*100
})

#longest inter-note number
lonIntNotNr = lapply(seg, function(x){
  x$syllable[x$pauseLen==max(x$pauseLen, na.rm = T)][1]
})

#shortest inter-note number
shoIntNotNr = lapply(seg, function(x){
  x$syllable[x$pauseLen==min(x$pauseLen, na.rm = T)][1]
})

###########################################################################################
#  save data as CSV file

## put all relevant variables in order into a data.frame, taking care to convert frequencies to Hz
df = cbind(recordID, calNr, "calStaTim" = rep("", length(fileName)), 
           "domFre"=unlist(domFre)*1000, "basFre"= unlist(basFre)*1000, basFreAmpRed, 
           "empFre2"= unlist(empFre2)*1000, empFre2AmpRed, 
           "banWidLow"=unlist(banWidLow)*1000,"banWidUpp"=unlist(banWidUpp)*1000,
           "domFreNotMd" =unlist(domFreNotMd)*1000,"domFreNotRan" =unlist(domFreNotRan)*1000,
           "domFre1Dev"=unlist(domFre1Dev)*1000,
           nNot, calDur, notRat, notDurMd, notDurRan, notDur1Dev, lonNotNr, shoNotNr,
           intNotIntMd, intNotIntRan, lonIntNotNr, shoIntNotNr)


###########################################################################################
#  extra information

## list of field names
fieldNames = colnames(df)

## information on field names
fieldInformation = c("recording ID", "call number", "call start time (min)", 
                     "dominant frequency (Hz)", "basal frequency (Hz)",
                     "basal frequenct ampluitude reduction (% of domFre)",
                     "2nd emphasised frequency (Hz)",
                     "2nd emph freq amplitude reduction (% of domFre)",
                     "lower 90% bandwidth peak (Hz)",
                     "upper 90% bandwidth peak (Hz)",
                     "dominant frequency median among notes (Hz)",
                     "range (max-min) of dominant frequency (Hz) among notes",
                     "1st note deviance in dominant frequency (Hz)",
                     "number of notes", "call duration (s)",
                     "note rate (notes/s)", 
                     "median note duration (ms)",
                     "range in note duration (% of median)",
                     "1st note deviance in duration (% of median)",
                     "longest note nr",
                     "shortest note nr",
                     "inter-note interval median (ms)",
                     "inter-note interval range (% of median)",
                     "longest inter-note interval nr",
                     "shortest inter-note interval nr")

df2 = cbind(fieldNames,fieldInformation)

require(qpcR)

## BIND THE two data frames, filling blanks with NA
df3 = qpcR:::cbind.na(df,df2)

## replace NA with ""
df3[is.na(df3)] = ""

##save as CSV 
write.csv(df3, paste0(getwd(),"/metadata per call.csv"), row.names = F)

#################################################################################
##write additional per note metadata for later use

## metadata per note
write.csv(cbind("recCalID" = unlist(lapply(1:length(seg), function(x){
  rep(calID[[x]], nrow(seg[[x]]))
})), 
"noteNr" = unlist(lapply(seg, function(x){
  1:nrow(x)
})),
"domFreNot" = unlist(domFreNot), 
"basFreNot" = unlist(basFreNot), 
"empFre2Not" = unlist(empFre2Not),
"notDur" = unlist(lapply(seg, function(x){
  x$sylLen
})), 
"intNotDur" = unlist(lapply(seg, function(x){
  x$pauseLen
}))),
paste0(getwd(),"/metadata per note.csv"), row.names = F, quote = F)

######################################################################################
## clear the memory
rm(list = ls())