R/seewave.r

Defines functions stdft sspectro soscillo sumsmooth rescale outputw inputw gaussian.w cunwrap zcr zc zapsilw write.audacity wf wav2flac wasp TKEO TFSD timer timelapse th synth2 aux symba spectro3D spectro specflux spec soundscapespec songmeterdiag songmeter smoothw simspec sh sfm setenv seedata sddB savewav rugo roughness rms rmnoise rmoffset rmam resamp revw repw read.audacity pulsew preemphasis playlist phaseplot2 phaseplot oscilloST oscillo octaves notefreq noisew NDSI mutew moredB micsens melfilterbank mel meanspec meandB M lts logspec.dist localpeaks lfs listen ks.dist kl.dist itakura.dist istft ifreq hilbert H ggspectro gammatone ftwindow fpeaks fma fir field ffilter fdoppler fbands fadew export env echo dynspectro dynspec dynoscillo duration drawfilter drawenv diffwave diffspec diffenv diffcumspec dfreq deletew dBweight cutw cutspec csh crest covspectro corspec corenv convSPL combfilter coh cepstro ceps ccoh bwfilter beep autoc audiomoth attenuation AR audiomoth.rename audiomoth ama akamatsu afilter addsilw

Documented in addsilw afilter akamatsu ama AR attenuation audiomoth audiomoth.rename autoc aux beep bwfilter ccoh ceps cepstro coh combfilter convSPL corenv corspec covspectro crest csh cunwrap cutspec cutw dBweight deletew dfreq diffcumspec diffenv diffspec diffwave drawenv drawfilter duration dynoscillo dynspec dynspectro echo env export fadew fbands fdoppler ffilter field fir fma fpeaks ftwindow gammatone gaussian.w ggspectro H hilbert ifreq inputw istft itakura.dist kl.dist ks.dist lfs listen localpeaks logspec.dist lts M meandB meanspec mel melfilterbank micsens moredB mutew NDSI noisew notefreq octaves oscillo oscilloST outputw phaseplot phaseplot2 playlist preemphasis pulsew read.audacity repw resamp rescale revw rmam rmnoise rmoffset rms roughness rugo savewav sddB seedata setenv sfm sh simspec smoothw songmeter songmeterdiag soscillo soundscapespec spec specflux spectro spectro3D sspectro stdft sumsmooth symba synth2 TFSD th timelapse timer TKEO wasp wav2flac wf write.audacity zapsilw zc zcr

################################################################################
##  seewave by Jerome Sueur, Caroline Simonis & Thierry Aubin
##  http://rug.mnhn.fr/seewave/
################################################################################

################################################################################
## ACI
################################################################################

ACI <- function (wave, f, channel = 1, wl = 512, ovlp = 0, wn = "hamming", flim = NULL, nbwindows = 1) 
	{
		input <- inputw(wave = wave, f = f, channel = channel)
		wave <- input$w
		f <- input$f
		rm(input)
		z <- sspectro(wave, f, wl = wl, wn = wn, ovlp = ovlp)
		if (!is.null(flim))
		{
			flim <- flim * 1000 * wl/f
			z <- z[flim[1]:flim[2], ]
		}
		
		acis <- array(0, nbwindows)
		
		for (j in 1:nbwindows) 
		{
			l <- dim(z)[2]
			z1 <- z[, floor(l/nbwindows * (j - 1) + 1):floor(l/nbwindows * j)]
			t <- array(0, dim(z1) - c(0, 1)) 
				for (i in 1:(dim(t)[1]))
				{ 
					mp <- (diff(z1[i, ]))/sum(z1[i, ])
                                        if(any(mp=='-1') | any(mp=='NaN')){t[i, ] <- as.numeric(NaN)}
                                        else {t[i, ] <- abs(mp)}
				}
			acis[j] <- sum(t)
		}
		return(sum(na.omit(as.vector(acis))))
	}



################################################################################
##                                ACOUSTAT
################################################################################

acoustat <- function (
            wave,
            f,
            channel = 1,
            wl = 512,
            ovlp = 0,
            wn = "hanning",
            tlim = NULL, 
            flim = NULL, 
            aggregate = sum,
            fraction = 90,
            plot = TRUE,
            type = "l",
            ...) 

{
    ## INPUT
    input <- inputw(wave = wave, f = f, channel = channel)
    wave <- input$w
    f <- input$f
    if (!is.null(tlim)) {
        wave <- cutw(wave, f = f, from = tlim[1], to = tlim[2])
    }
    n <- length(wave)
    rm(input)

    ## STFT
    m <- sspectro(wave, f = f, wl = wl, ovlp = ovlp, wn = wn)

    ## FREQUENCY SELECTION AND FREQUENCY AXIS
    if (is.null(flim)) {freq <- seq(0, (f/2) - (f/wl), length.out = nrow(m))/1000}
    else {
        fl1 <- flim[1] * nrow(m) * 2000 / f
        fl2 <- flim[2] * nrow(m) * 2000 / f
        m <- m[(fl1:fl2) + 1, ]
        freq <- seq(flim[1], flim[2], length.out = nrow(m))
    }

    ## TIME AXIS
    time <- seq(0, n/f, length.out = ncol(m))

    ## CONTOURS
    t.cont <- apply(m, MARGIN = 2, FUN = aggregate)
    f.cont <- apply(m, MARGIN = 1, FUN = aggregate)
    t.cont <- t.cont/sum(t.cont)
    f.cont <- f.cont/sum(f.cont)
    t.cont.cum <- cumsum(t.cont)
    f.cont.cum <- cumsum(f.cont)

    ## STATISTICS
    P <- fraction/100
    proportions <- as.matrix(c((1 - P)/2, 0.5, P + (1 - P)/2))
    t.quantiles <- apply(proportions, MARGIN = 1, function(x) time[length(t.cont.cum[t.cont.cum <= 
        x]) + 1])
    f.quantiles <- apply(proportions, MARGIN = 1, function(x) freq[length(f.cont.cum[f.cont.cum <= 
        x]) + 1])
    time.P1 <- t.quantiles[1]
    time.M <- t.quantiles[2]
    time.P2 <- t.quantiles[3]
    time.IPR <- time.P2 - time.P1
    freq.P1 <- f.quantiles[1]
    freq.M <- f.quantiles[2]
    freq.P2 <- f.quantiles[3]
    freq.IPR <- freq.P2 - freq.P1
    results <- list(time.contour = cbind(time = time, contour = t.cont), 
        freq.contour = cbind(frequency = freq, contour = f.cont), 
        time.P1 = time.P1, time.M = time.M, time.P2 = time.P2, 
        time.IPR = time.IPR, freq.P1 = freq.P1, freq.M = freq.M, 
        freq.P2 = freq.P2, freq.IPR = freq.IPR)

    ## PLOT AND RETURN
    if (plot) {
        def.par <- par(no.readonly = TRUE)
        on.exit(par(def.par))
        par(mfrow = c(2, 1))
        plot(x = time, y = t.cont, type = type, xlab = "Time (s)", 
            ylab = "Probability", main = "Time envelope", ...)
        segments(x0 = t.quantiles[1], y0 = 0, y1 = t.cont[time == 
            time.P1], col = 2)
        segments(x0 = t.quantiles[2], y0 = 0, y1 = t.cont[time == 
            time.M], col = 2, lwd = 2)
        segments(x0 = t.quantiles[3], y0 = 0, y1 = t.cont[time == 
            time.P2], col = 2)
        plot(x = freq, y = f.cont, type = type, xlab = "Frequency (kHz)", 
            ylab = "Probability", main = "Frequency spectrum", 
            ...)
        segments(x0 = f.quantiles[1], y0 = 0, y1 = f.cont[freq == 
            freq.P1], col = 2)
        segments(x0 = f.quantiles[2], y0 = 0, y1 = f.cont[freq == 
            freq.M], col = 2, lwd = 2)
        segments(x0 = f.quantiles[3], y0 = 0, y1 = f.cont[freq == 
            freq.P2], col = 2)
        invisible(results)
    }
    else return(results)
}



################################################################################
##                                ADDSILW
################################################################################

addsilw<-function(
                  wave,
                  f,
                  channel = 1,
                  at = "end",
                  choose = FALSE,
                  d = NULL,
                  plot = FALSE,
                  output = "matrix",
                  ...)

{
  input<-inputw(wave=wave,f=f, channel=channel); wave<-input$w ; f<-input$f ; bit <- input$bit ; rm(input)

  switch(at,
         start = {at<-0},
         middle = {at<-nrow(wave)/(2*f)},
         end = {at<-nrow(wave)/f}
         )
  
  if(is.null(d)) stop("silence duration has to be set with the argument 'd'")

  if(choose)
    { 
      cat("choose position on the wave\n")
      if(.Platform$OS.type == "windows") flush.console()
      oscillo(wave,f=f)
      coord<-locator(n=1)
      at<-coord$x[1]; abline(v=at,col=2,lty=2)
    }

  pos <- round(at*f)
  wave1 <- wave[1:pos,1]
  sil <- rep(0,d*f)
  wave2 <- wave[-(1:pos),1]

  wave3 <- c(wave1,sil,wave2)

  wave3 <- outputw(wave=wave3, f=f, format=output)
  
  if(plot)
    {
      oscillo(wave=wave3,f=f,...)
      invisible(wave3)
    }
  else return(wave3)
}

################################################################################
##                                AFILTER                                       
################################################################################

afilter<-function(
                  wave,
                  f,
                  channel = 1,
                  threshold = 5,
                  plot = TRUE,
                  listen = FALSE,                  
                  output = "matrix",
                  ...
                  )

{
  input<-inputw(wave=wave,f=f,channel=channel) ; wave<-input$w ; f<-input$f ; bit <- input$bit ; rm(input)

  t1<-max(abs(wave))*(threshold/100)
  wave1 <- ifelse(abs(wave)<=t1,yes=0,no=1)
  wave2 <- wave[,1]*wave1[,1]
  wave2 <- outputw(wave=wave2, f=f, format=output)

  if(plot)
    {
      oscillo(wave=wave2,f=f,...)
      invisible(wave2)
    }
  else
    {
      if(listen) {listen(wave2,f=f)}
      return(wave2)
    }
}


################################################################################
##                                AKAMATSU
################################################################################


akamatsu <- function(Lx,         # L in cm
                     Ly,
                     Lz,
                     mode = c(1,1,1),
                     c = 148000,    # in cm/s
                     plot = FALSE,
                     xlab = "Frequency (kHz)",
                     ylab = "Attenuation distance (cm)",
                     ...)

{          
## if loop to be sure that the largest length is the right one
    if (Lx > Ly){
        a  <- Lx
        Lx <- Ly
        Ly <- a
    }

    Fres <- (c*sqrt((mode[1]/Lx)^2+(mode[2]/Ly)^2+(mode[3]/Lz)^2))/2 # resonance frequency of the aquarium 
    Fcut <- (c*sqrt((1/Ly^2)+(1/Lz^2)))/2 # frequency of the high pass filter 
    F <- list(res=Fres, cut=Fcut)

## Plot showing the attenuation distance in function of the frequency 
    if(plot == TRUE){
        f <- seq(0, Fcut, Fcut/100)
        f <- f[-length(f)]
        D <- (2*log(10)*c)/(4*pi*Fcut*sqrt(1-(f/Fcut)^2))
        plot(f/1000, D, xlab = xlab, ylab = ylab,...)
        invisible(F)
    }

    return(F)
}


################################################################################
##                                AMA
################################################################################

ama<-function(
              wave,
              f,
              channel = 1,
              envt = "hil",
              wl = 512,
              plot = TRUE,
              type = "l",
              ...)

{
  input<-inputw(wave=wave,f=f,channel=channel) ; wave<-input$w ; f<-input$f ; rm(input)
  enve<-env(wave=wave, f=f, envt = envt, plot = FALSE)
  meanspec(wave=enve, f=f, wl=wl, plot=plot, type=type,...)
}


################################################################################
##                                AUDIOMOTH
################################################################################

audiomoth <- function(x,                 # a character vector, not a Wave object
                      tz = ""            # a character vector defining a time zone specification, see as.POSIXct(), argument 'tz'
                      )
  {
    ## INPUT
    if (!is.character(x)) stop("'x' should be of mode character.")
 
    ## PREPARE RESULTS
    options(stringsAsFactors = FALSE)
    res <- data.frame(year = numeric(0), month = numeric(0), day = numeric(0),
                      hour = numeric(0), min = numeric(0), sec = numeric(0),
                      time = numeric(0))
    N <- length(x)  ## number of file names    
    ## LOOP
    for (i in 1:N)
    {
      tmp <- x[i]
      n <- nchar(tmp)
      extension <- substr(tmp, start = n-3, stop = n)
      if(extension != ".wav" & extension != ".WAV"){warning(paste("File '", tmp, "' is not a '.wav' file", sep=""))}
      else
      {
        hex <- unlist(strsplit(tmp, extension))
        num <- strtoi(hex, base = 16)
        time <- as.POSIXct(num, tz=tz, origin = "1970-01-01")  ## POSIXct as in songmeter()
        year <- as.numeric(format(time, "%Y"))
        month <- as.numeric(format(time, "%m"))
        day <- as.numeric(format(time, "%d"))
        hour <- as.numeric(format(time, "%H"))
        min <- as.numeric(format(time, "%M"))
        sec <- as.numeric(format(time, "%S"))
        res <- rbind(res, data.frame(year, month, day, hour, min, sec, time))
      }
    }
    return(res)    
    options(stringsAsFactors = TRUE)
}


################################################################################
##                                AUDIOMOTH.RENAME
################################################################################

audiomoth.rename <- function(
                             dir,
                             overwrite = FALSE,
                             tz = "",
                             prefix = ""
                             )
{
    if(!is.character(dir)) stop("'dir' should be a character vector giving the path to a directory")
    files <- dir(dir, pattern="[WAV]$")
    from <- audiomoth(files, tz=tz)
    if(prefix!="") prefix  <- paste(prefix, "_", sep="")
    out <- function(x) formatC(x, flag="0", digits=1, format="d")
    to <- paste(prefix, out(from$year), out(from$month), out(from$day), "_", out(from$hour), out(from$min), out(from$sec), ".wav", sep="")
    from <- paste(dir, files, sep="/")
    to <- paste(dir, to, sep="/")
    if(overwrite) {res <- file.rename(from=from, to=to)} else {res <- file.copy(from=from, to=to)}
    return(res)
}

################################################################################
## AR                                        
################################################################################

AR <- function(...,
               datatype = "objects",    # if "file" -> only .wav and .mp3 read by tuneR functions readWave and readMP3
               envt = "hil",
               msmooth = NULL,
               ksmooth = NULL,
               ssmooth = NULL,
               pattern = "[wav]$|[WAV]$|[mp3]$")
    
{
    ## data as R objects
    if(datatype=="objects"){
        if(is.character(list(...)[[1]])) stop("The first argument should not be a character string if 'datatype=object'")
        ## number of objects = number of arguments - missing arguments
        n <- nargs() - !missing(datatype) - !missing(envt) - !missing(pattern) 
        ## number of objects = number of arguments - NULL arguments
        if(is.null(msmooth)) m <- 0 else m <- 1
        if(is.null(ksmooth)) k <- 0 else k <- 1
        if(is.null(ssmooth)) s <- 0 else s <- 1
        n <- n - m - k- s
        if(n <= 1) warning("Computing the AR index on a single object does not make a lot of sense as AR is a ranked index")
        name <- M <- Ht <- AR <- numeric(n)
        for(i in 1:n){ # does not use the function M() to avoid to compute the envelope twice
            wave <- list(...)[[i]]
            if(class(wave)[1]=="Wave"|class(wave)[1]=="WaveMC") {bit <- wave@bit}
            else if(class(wave)[1]=="audioSample") {bit <- wave$bits}
            if(!is.null(names(list(...))[i])) {name[i] <- names(list(...))[i]}
            env <- env(wave, envt=envt, msmooth=msmooth, ksmooth=ksmooth, ssmooth=ssmooth, plot=FALSE)
            env.norm <- env/sum(env)
            M[i] <- median(env)*2^(1-bit)
            Ht[i] <- th(env.norm)
        }
    }
    ## data as .wav or .mp3 files
    ## data should be a path to a directory containing either .wav or .mp3 files
    else if(datatype=="files"){
        if(!is.character(...)) stop("The first argument should be a character string when 'datatype='files''")
        files <- dir(..., pattern=pattern)
        n <- length(files)
        if(n <= 1) warning("Computing the AR index on a single object does not make a lot of sense as AR is a ranked index")
        name <- M <- Ht <- AR <- numeric(n)
        for(i in 1:n){ # does not use the function M() to avoid to compute the envelope twice
            filename <- basename(files[i])
            nch <- nchar(filename)
            extension <- substr(filename, start = nch-3, stop = nch)
            name[i] <- substr(filename, start=1, stop=nch-4)
            if(extension == ".wav" | extension == ".WAV") {wave <- readWave(files[i])}
            else if (extension == ".mp3") {wave <- readMP3(files[i])}
            env <- env(wave, envt=envt, msmooth=msmooth, ksmooth=ksmooth, ssmooth=ssmooth, plot=FALSE)
            env.norm <- env/sum(env)
            M[i] <- median(env)*2^(1-wave@bit)
            Ht[i] <- th(env.norm)
        }
    }
    
    ## OUTPUT
    AR <- (rank(M)*rank(Ht))/(n^2)
    res <- data.frame(M, Ht, AR)
    if(!is.numeric(name)) rownames(res) <- name
    return(res)
}


################################################################################
##                                ATTENUATION                                        
################################################################################

attenuation <- function(lref,
                        dref = 1,
                        dstop,
                        n,
                        plot = TRUE,
                        xlab = "Distance (m)",
                        ylab = "dB",
                        type="l",
                        ...
                        )

{
    d <- seq(from=dref, to=dstop, length.out=n)
    dB <- lref-(20*log10(d/dref))

    if(plot)
        {
            plot(x=d, y=dB, xlab=xlab, ylab=ylab, type=type,...)
            invisible(dB)
        }
    else {return(dB)}
}


################################################################################
## AUDIOMOTH
################################################################################

audiomoth <- function(x,                 # a character vector, not a Wave object
                      tz = ""            # a character vector defining a time zone specification, see as.POSIXct(), argument 'tz'
                      )
  {
    ## INPUT
    if (!is.character(x)) stop("'x' should be of mode character.")
 
    ## PREPARE RESULTS
    options(stringsAsFactors = FALSE)
    res <- data.frame(year = numeric(0), month = numeric(0), day = numeric(0),
                      hour = numeric(0), min = numeric(0), sec = numeric(0),
                      time = numeric(0))
    N <- length(x)  ## number of file names    
    ## LOOP
    for (i in 1:N)
    {
      tmp <- x[i]
      n <- nchar(tmp)
      extension <- substr(tmp, start = n-3, stop = n)
      if(extension != ".wav" & extension != ".WAV") {warning(paste("File '", tmp, "' is not a '.wav' file", sep=""))}
      else
      {
        hex <- unlist(strsplit(tmp, extension))
        num <- strtoi(hex, base = 16)
        time <- as.POSIXct(num, tz=tz, origin = "1970-01-01")  ## POSIXct as in songmeter()
        year <- as.numeric(format(time, "%Y"))
        month <- as.numeric(format(time, "%m"))
        day <- as.numeric(format(time, "%d"))
        hour <- as.numeric(format(time, "%H"))
        min <- as.numeric(format(time, "%M"))
        sec <- as.numeric(format(time, "%S"))
        res <- rbind(res, data.frame(year, month, day, hour, min, sec, time))
      }
    }
    return(res)    
    options(stringsAsFactors = TRUE)
}


################################################################################
## AUTOC                                         
################################################################################

autoc <- function(
                wave,
                f,
                channel = 1,
                wl = 512,
                fmin = 1,
                fmax = f/2,
                threshold = NULL,
                plot = TRUE,
                xlab = "Time (s)",
                ylab = "Frequency (kHz)",
                ylim = c(0,f/2000),
                pb = FALSE,
                ...)

{
  input <- inputw(wave=wave,f=f,channel=channel) ; wave<-input$w ; f<-input$f ; rm(input)

  if(!is.null(threshold)) wave <- afilter(wave=wave,f=f,threshold=threshold,plot=FALSE)

  lag.min <- round(wl*(fmin/(f/2)))
  lag.max <- round(wl*(fmax/(f/2)))
  if(lag.max==wl) {lag.max <- lag.max-1}
  n <- nrow(wave)
  step <- seq(1,n-2*wl,wl)
  N <- length(step) 
  if(pb) {pbar <- txtProgressBar(min=0, max=N, style = 3)}
  R <- matrix(data=numeric((lag.max+1)*N), lag.max+1, N)
  for (i in step)
    {
      R[,which(step==i)] <- as.vector(acf(wave[i:(wl+i-1)], lag.max=lag.max, plot=FALSE)$acf)
    }
  R[is.nan(R)|is.na(R)] <- 0 # acf() can produce NaN or NA
  tfond <- numeric(N)
  excl <- 1:lag.min

  for (i in 1:N)
    {
    tfond[i] <- fpeaks(R[-excl,i], f=NA, nmax=1)[1]
    if(pb) {setTxtProgressBar(pbar, i)}
    }

  tfond <- tfond + length(excl)
  y<-f/tfond/1000
  x<-seq(0,n/f,length.out=N)  

  if(plot)
    {
      plot(x=x, y=y,
           xlab = xlab,
           ylab = ylab, ylim = ylim,
           las = 1,
           ...)
      invisible(cbind(x,y))
    }
  else {return(cbind(x,y))}
}


################################################################################
##  BEEP                                       
################################################################################

beep <- function(d=0.5, f=8000, cf=1000){
    s <- synth(d=d, f=f, cf=cf, output="matrix")
    listen(s, f=f)
}


################################################################################
##  BWFILTER                                       
################################################################################

bwfilter <- function(
                    wave,
                    f,
                    channel = 1,
                    n = 1,
                    from = NULL,
                    to = NULL,
                    bandpass = TRUE,
                    listen = FALSE,
                    output = "matrix"
                    )
    
{
  ## input
  input<-inputw(wave=wave,f=f,channel=channel) ; wave<-input$w ; f<-input$f ; bit <- input$bit ; rm(input)

  ## Butterworth filter
  if(isTRUE(bandpass)){
    if(is.null(from)  && !is.null(to)) type <- "low"  ## low-pass
    if(!is.null(from) && is.null(to))  type <- "high" ## high-pass
    if(!is.null(from) && !is.null(to)) type <- "pass" ## bandpass
  }
  else{
    if(is.null(from)  && !is.null(to)) type <- "high" ## high-pass 
    if(!is.null(from) && is.null(to))  type <- "low"  ## low-pass
    if(!is.null(from) && !is.null(to)) type <- "stop" ## bandstop
  }
  
  bf <- signal::butter(n=n, W=c(from,to)/(f/2), type=type)
  
  ## filtration
  wave2 <- signal::filtfilt(filt=bf, x=wave)
  
  ## output
  wave2 <- outputw(wave=wave2, f=f, format=output)
  if(listen) {listen(wave2,f=f)}
  return(wave2)
}


################################################################################
##                                CCOH                                        
################################################################################

ccoh<-function(
               wave1,
               wave2,
               f,
               channel = c(1,1),
               wl = 512,
               ovlp = 0,
               plot = TRUE,
               grid = TRUE,
               scale = TRUE,
               cont = FALSE,
               collevels = seq(0,1,0.01),
               palette = reverse.heat.colors,
               contlevels = seq (0,1,0.01),
               colcont = "black",
               colbg = "white",
               colgrid = "black",
               colaxis = "black",
               collab = "black",
               xlab = "Time (s)",
               ylab = "Frequency (kHz)",
               scalelab = "Coherence",
               main = NULL, 
               scalefontlab = 1,
               scalecexlab =0.75,
               axisX = TRUE,
               axisY = TRUE,
               flim = NULL,
               flimd = NULL,
               ...)

{
  if(length(channel)!=2) stop("'channel' should be a numeric vector of length 2")
  input1<-inputw(wave=wave1,f=f,channel=channel[1]) ; wave1<-input1$w ; f<-input1$f ; rm(input1)
  wave2<-inputw(wave=wave2,f=f,channel=channel[2])$w

  n1<-nrow(wave1)
  n2<-nrow(wave2)
  if(n1 != n2) stop("'wave 1' and 'wave 2' must have the same length")
  n<-n1


                                        # dynamic vertical zoom (modifications of analysis parameters)
  if(!is.null(flimd))
    {
                                        # zoom magnification
      mag<-round((f/2000)/(flimd[2]-flimd[1]))
                                        # new parameters
      wl<-wl*mag
      if(ovlp==0) ovlp<-100
      ovlp<-100-round(ovlp/mag)
                                        # use of normal flim for following axis modifications
      flim<-flimd
    }

  step<-seq(1,n+1-wl,wl-(ovlp*wl/100))	# coherence windows, +1 added @ 2017-04-20

  z1<-matrix(data=numeric((wl)*length(step)),wl,length(step))

  for(i in step)
    {
      z1[,which(step==i)]<-spec.pgram(cbind(wave1[i:(wl+i-1),],
                  wave2[i:(wl+i-1),]), spans = c(3,3), fast=FALSE, taper=FALSE, plot=FALSE)$coh
    }

  z<-z1[1:(wl/2),]							  

                                        # X axis settings
  X<-seq(0,n/f,length.out=length(step))

                                        # vertical zoom
  if(is.null(flim))
    {
      Y<-seq(0,f/2000,length.out=nrow(z))
    }
  else
    {
      fl1<-flim[1]*nrow(z)*2000/f
      fl2<-flim[2]*nrow(z)*2000/f
      z<-z[fl1:fl2,]
      Y<-seq(flim[1],flim[2],length.out=nrow(z))
    }
  
  Z<-t(z)
  
  if(plot)
    {
      Zlim<-range(Z, finite = TRUE) 
      
      if(scale)
        {
          def.par <- par(no.readonly = TRUE)
          on.exit(par(def.par))
          layout(matrix(c(1, 2), ncol = 2, byrow=TRUE), widths = c(6, 1))
          par(mar=c(5,4.1,1,0),las=1,cex=1,bg=colbg,col=colaxis,col.axis=colaxis,col.lab=collab)
          filled.contour.modif2(x=X ,y=Y, z=Z, levels=collevels, nlevels=20,
                                plot.title=title(main=main,xlab=xlab,ylab=ylab), color.palette=palette,axisX=axisX, axisY=axisY)
          if(grid) grid(nx=NA, ny=NULL, col=colgrid)
          if(colaxis != colgrid) abline(h=0,col=colaxis) else abline(h=0,col=colgrid)
          par(mar=c(5,1,4.5,3),las=0)
          dBscale(collevels=collevels,palette=palette,fontlab=scalefontlab,
                  cexlab=scalecexlab,collab=collab,textlab=scalelab,colaxis=colaxis)
        }
      
      if(scale==FALSE)
        {
          par(las=1, col=colaxis, col.axis=colaxis, col.lab=collab,,bg=colbg,...)
          filled.contour.modif2(x=X ,y=Y, z=Z, levels=collevels, nlevels=20,
                                plot.title=title(main=main,xlab=xlab,ylab=ylab), color.palette=palette, axisX=axisX, axisY=axisY,
                                col.lab=collab,colaxis=colaxis)		
          if(grid) grid(nx=NA, ny=NULL, col=colgrid)
          if(colaxis != colgrid) abline(h=0,col=colaxis) else abline(h=0,col=colgrid)
        }

      if(cont) 
        {
          contour(X,Y,Z,add=TRUE,
                  levels=contlevels,nlevels=5,col=colcont,...)
        }  
      invisible(list(time=X, freq=Y, coh=Z))
    }
  else return(list(time=X, freq=Y, coh=Z))
}




################################################################################
##                                CEPS                                         
################################################################################

ceps <- function(
                 wave,
                 f,
                 channel = 1,
                 phase = FALSE,
                 wl = 512,
                 at = NULL,
                 from = NULL,
                 to = NULL,
                 tidentify = FALSE,   # identify in seconds
                 fidentify = FALSE,   # identify in Hz
                 col = "black",
                 cex = 1,
                 plot = TRUE,
                 qlab = "Quefrency (bottom: s, up: Hz)",
                 alab = "Amplitude",
                 qlim = NULL,
                 alim = NULL,
                 type= "l",
                 ...
                 )

{
  input<-inputw(wave=wave,f=f,channel=channel) ; wave<-input$w ; f<-input$f ; rm(input)

  if(!is.null(from)|!is.null(to))
    {
      if(is.null(from) && !is.null(to)) {a<-1; b<-round(to*f)}
      if(!is.null(from) && is.null(to)) {a<-round(from*f); b<-length(wave)}
      if(!is.null(from) && !is.null(to))
        {
          if(from>to) stop("'from' cannot be superior to 'to'")
          if(from==0) {a<-1} else a<-round(from*f)
          b<-round(to*f)
        }
      wave<-as.matrix(wave[a:b,])
    }

  if(!is.null(at))
    {
      if(wl==FALSE) stop("Argument 'wl' has to be set up, for instance wl=512")
      c<-round(at*f)
      wl2<-wl%/%2
      wave<-as.matrix(wave[(c-wl2):(c+wl2),])
    }

  n<-nrow(wave)
  N<-round(n/2)

  h <- fft(wave[,1])                                                    
  logh <- log(abs(h))
  if(phase) {logh <- logh + complex(real=0, imaginary=1)*cunwrap(Arg(h))} 
  z <- Re(fft(logh, inverse=TRUE))                                  
  z <- z[1:N]                                           
  x <- seq(0,N/f,length.out=N)                           

  if(plot)
    {
      plot(x=x,y=z,
           xlab=qlab,xaxt="n",xaxs="i",
           ylab=alab,yaxt="n",yaxs="i",
           type=type,col=col,cex=cex,xlim=qlim,ylim=alim,...)
      if(!is.null(qlim)) E<-qlim[2] else E<-N/f
      X<-seq(0,E,length.out=7)
      axis(side=1,at=X, labels=round(X,3))
      axis(side=3,at=X, labels=round(1/X,3))

      if(tidentify)
        {
          cat("time identification: choose points on the cepstrum\n")
          if(.Platform$OS.type == "windows") flush.console()
          id<-identify(x=x,y=z,labels=round(x,5),tolerance=0.15,col="red")
          return(round(x[id],5))
        }
      
      if(fidentify)
        {
          cat("frequency identification: choose points on the cepstrum\n")
          if(.Platform$OS.type == "windows") flush.console()
          id<-identify(x=x,y=z,labels=round(1/round(x,5),1),tolerance=0.15,col="red")
          return(round(1/round(x[id],5),1))
        }
      results <- cbind(x,z)
      invisible(results)
    }

  else
    {
      results <- cbind(x,z)
      return(results)
    }
}



################################################################################
##                                   CEPSTRO                                    
################################################################################

cepstro <- function(
                  wave,
                  f,
                  channel = 1,
                  wl = 512,
                  ovlp = 0,
                  plot = TRUE,
                  grid = TRUE,
                  scale = TRUE,
                  cont = FALSE,
                  collevels = seq(0,1,0.01),
                  palette = reverse.heat.colors,
                  contlevels = seq (0,1,0.01),
                  colcont = "black",
                  colbg = "white",
                  colgrid = "black",
                  colaxis = "black",
                  collab = "black",
                  xlab = "Time (s)",
                  ylab = "Quefrency (ms)",
                  scalelab = "Amplitude",
                  main = NULL,
                  scalefontlab = 1,
                  scalecexlab =0.75,
                  axisX = TRUE,
                  axisY = TRUE,
                  tlim = NULL,  
                  qlim = NULL, 
                  ...)

{
  input <- inputw(wave=wave,f=f,channel=channel) ; wave <- input$w ; f <- input$f ; rm(input)
  wave <- ifelse(wave==0,yes=1e-6,no=wave)

  if(!is.null(tlim)) wave <- cutw(wave,f=f,from=tlim[1],to=tlim[2])                                      

  n <- nrow(wave)
  p <- round(n/2)
  step <- seq(1,n+1-wl,wl-(ovlp*wl/100)) # +1 added @ 2017-04-20

  N <- length(step)
  WL <- wl%/%2
  z <- matrix(data=numeric(wl*N),wl,N)
  for(i in step)
    {z[,which(step==i)]<-Re(fft(log(abs(fft(wave[i:(wl+i-1),]))),inverse=TRUE))}
  z <- z[1:WL,]
  z <-ifelse(z=="NaN"|z=="-Inf"|z<=0|z=="Inf",yes=0,no=z)

 ## Time X-axis settings
  X <- seq(0,n/f,length.out=length(step))

 ## Quefrency Y-axis settings
  Y<-seq(0, WL/f, length.out=nrow(z))*1000
  if(!is.null(qlim)) {
      ql1 <- which(round(Y,1)==qlim[1])[1]
      ql2 <- which(round(Y,1)==qlim[2])[1]
      z <- z[ql1:ql2,]
      Y <- seq(qlim[1], qlim[2], length.out=nrow(z))
  }

## Amplitude Z-axis settings
  Z<-t(z/max(z))

  if(plot)
    {
      Zlim<-range(Z, finite = TRUE)
      if(scale)
        {
          layout(matrix(c(1, 2), ncol = 2, byrow=TRUE), widths = c(6, 1))
          par(mar=c(5,4.1,1,0),las=1,cex=1,col=colaxis,col.axis=colaxis,col.lab=collab,bg=colbg)
          filled.contour.modif2(x=X ,y=Y, z=Z,
                                levels=collevels, nlevels=20,
                                plot.title=title(main=main,xlab=xlab,ylab=ylab),
                                color.palette=palette,axisX=axisX, axisY=axisY)
          if(grid) grid(nx=NA, ny=NULL, col=colgrid)
          if(colaxis != colgrid) abline(h=0,col=colaxis) else abline(h=0,col=colgrid)
          par(mar=c(5,1,4.5,3),las=0)
          dBscale(collevels=collevels,palette=palette,fontlab=scalefontlab,
                  cexlab=scalecexlab,collab=collab,textlab=scalelab,colaxis=colaxis)
        }

      if(scale==FALSE)
        {
          par(las=1, col=colaxis, col.axis=colaxis, col.lab=collab,bg=colbg,...)
          filled.contour.modif2(x=X ,y=Y, z=Z, levels=collevels, nlevels=20,
                                plot.title=title(main=main,xlab=xlab,ylab=ylab),
                                color.palette=palette, axisX=axisX, axisY=axisY,
                                col.lab=collab,colaxis=colaxis)
          if(grid) grid(nx=NA, ny=NULL, col=colgrid)
          if(colaxis != colgrid) abline(h=0,col=colaxis) else abline(h=0,col=colgrid)
        }

      if(cont)
        {
          contour(X,Y,Z,add=TRUE,
                  levels=contlevels,nlevels=5,col=colcont,...)
        }
      invisible(list(time=X, quef=Y, amp=t(Z)))
    }

  else{return(list(time=X, quef=Y, amp=t(Z)))}
}



################################################################################
##                                COH                                         
###############################################################################

coh<-function(
              wave1,
              wave2,
              f,
              channel = c(1,1),
              plot =TRUE,
              xlab = "Frequency (kHz)",
              ylab = "Coherence",
              xlim = c(0,f/2000),
              type = "l",
              ...
              )

{
  if(length(channel)!=2) stop("'channel' should be a numeric vector of length 2")
  input1<-inputw(wave=wave1,f=f,channel=channel[1]) ; wave1<-input1$w ; f<-input1$f ; rm(input1)
  wave2<-inputw(wave=wave2,f=f,channel=channel[2])$w

  n1<-nrow(wave1)
  n2<-nrow(wave2)
  if(n1 != n2) stop("'wave 1' and 'wave 2' must have the same length")

  Y<-spec.pgram(cbind(wave1, wave2), fast=FALSE, taper=FALSE,
                spans = c(3,3), plot=FALSE)$coh
  X<-seq(0,f/2000,length.out=nrow(Y))

  if(plot)
    {
      plot(x=X,y=Y,xlab=xlab,ylab=ylab,xlim=xlim,type=type,...)
      invisible(cbind(X,Y))
    }
  else return(cbind(X,Y))
}


################################################################################
## COMBFILTER                                         
###############################################################################

combfilter <- function(wave,
                       f,
                       channel = 1,
                       alpha,
                       K,
                       units = c("samples", "seconds"),
                       plot = FALSE,
                       output="matrix",
                       ...)
{
    ## INPUT
    input <- inputw(wave=wave,f=f,channel=channel) ; wave <- input$w ; f <- input$f ; rm(input)
    wave <- wave/max(wave)
    n <- length(wave)
    K <- switch(match.arg(units), samples = K, seconds = round(K*f))

    ## FILTER
    wave1 <- c(rep(0,K), wave[1:(n-K)]) + alpha*wave
    wave2 <- outputw(wave=wave1, f=f, format=output)

    ## FREQUENCY RESPONSE
    w <- seq(0, pi, length.out=n)
    H <- sqrt(1+alpha^2 + 2*alpha*cos(w*K))
    
    ## PLOT AND RETURN
    if(plot){
       def.par <- par(no.readonly = TRUE)
       on.exit(par(def.par))
       layout(matrix(c(1, 2), ncol = 2, byrow=TRUE), widths = c(4, 1))
        spectro(wave2, f=f, scale=FALSE, osc=FALSE, ...)
        par(mar=c(5.1,0,4.1,1))
        freq <- seq(0, f/2000, length.out=n)
        plot(x=H, y=freq, xlim=c(0, max(H)), type="l", xaxs="i", yaxs="i", xlab="Linear amplitude", ylab="", yaxt="n")
        invisible(wave2)
    }
    else return(wave2)
}


################################################################################
##                                CONVSPL                                         
###############################################################################

convSPL<-function(
                  x,
                  d = 1,
                  Iref = 10^-12,
                  pref = 2*10^-5 
                  )

{
  P<-4*pi*(d^2)*Iref*(10^(x/10))
  I<-Iref*(10^(x/10))
  p<-pref*(10^(x/20))
  conv<-list(P = P, I = I, p = p)
  return(conv)
}


################################################################################
##                                CORENV                                         
################################################################################

corenv <- function(
                   wave1,
                   wave2,
                   f,
                   channel = c(1,1),
                   envt="hil",
                   msmooth = NULL,
                   ksmooth = NULL,
                   ssmooth = NULL,
                   plot = TRUE,
                   plotval = TRUE,
                   method = "spearman",
                   col = "black",
                   colval = "red",
                   cexval = 1,
                   fontval = 1,
                   xlab = "Time (s)",
                   ylab = "Coefficient of correlation (r)",
                   type = "l",
                   pb = FALSE,
                   ...)

{
  options(warn=-1) ## to remove warning messages from cor.test()
  if(length(channel)!=2) stop("'channel' should be a numeric vector of length 2")
  input1<-inputw(wave=wave1,f=f,channel=channel[1]) ; wave1<-input1$w ; f<-input1$f ; rm(input1)
  wave2<-inputw(wave=wave2,f=f,channel=channel[2])$w

  n<-nrow(wave1)

#  cat("please wait...\n")
#  if(.Platform$OS.type == "windows") flush.console()
  
  x <- env(wave=wave1,f=f,envt=envt,msmooth=msmooth,ksmooth=ksmooth,ssmooth=ssmooth, plot=FALSE)
  y <- env(wave=wave2,f=f,envt=envt,msmooth=msmooth,ksmooth=ksmooth,ssmooth=ssmooth, plot=FALSE)

  nx <- nrow(x)
  ny <- nrow(y)

  if(nx != ny) stop("'wave 1' and 'wave 2' must have the same length")

  meanx <- mean(x)
  meany <- mean(y)
  diffx <- x-meanx
  r1 <- numeric(nx)
  r2 <- numeric(nx)

  if(pb) {pbar <- txtProgressBar(min=0, max=nx-1, style = 3)}

  for (i in 0:(nx-1))
    {
    r1[i+1] <- cor(x=x,y=c(y[(i+1):ny],rep(0,i)),method = method)
    r2[i+1] <- cor(x=x,y=c(rep(0,i),y[1:(ny-i)]),method = method)
    if(pb) {setTxtProgressBar(pbar, i)}
    }

  r2 <- r2[-1]
  r <- c(rev(r1),r2)
  rmax <- max(r,na.rm=TRUE)
  offset <- which.max(r)
  if(offset<=(length(r)/2)) {offsetp <- which.max(r)} else {offsetp <- which.max(r)-1}
  if(!is.null(msmooth)|!is.null(ksmooth)) {F <- f*nx/n; offsett <- (offsetp-nx)/F}
  else{offsett <- (offsetp-n)/f}

  if(offsetp < nx){p <- cor.test(x=x, y=c(y[(nx-offsetp+1):ny],rep(0,nx-offsetp)),method = method)}
  else { p<- cor.test(x=x, y=c(rep(0,offsetp-nx),y[1:(ny-(offsetp-nx))]), method = method)}
  p <- p$p.value

  X <- seq(-n/f,n/f,length.out=2*nx-1)
  corr <- list(r = cbind(X,r), rmax = rmax, p = p, t = offsett)

  options(warn=0) ## to reset warning messages

  if(plot)
    {
      plot(x = X, y = r, xlab = xlab, ylab = ylab, col = col, type = type,...)
      if(plotval)
        {
          mtext(paste(
                      "rmax = ", as.character(round(rmax,2)),
                      ", offset = ", as.character(round(offsett,3)), "s", sep=" "),
                side=3, line=-2, col=colval, cex=cexval, font=fontval)
          segments(x0=round(offsett,3), y0=min(r), x1=round(offsett,3), y1=rmax, col=colval, lty=2)
          segments(x0=-n/f, y0=rmax, x1=round(offsett,3), y1=rmax, col=colval, lty=2)
          points(x=round(offsett,3), y=rmax, pch=19, cex=1, col=colval)
  
        }
      return(invisible(corr))    # 2016-12-01 does not work if not wrapped by return
    }
  else
    {
      return(corr)
    }
 if(pb) close(pbar)
}



################################################################################
##                                CORSPEC                                         
################################################################################

corspec <- function(
                    spec1,
                    spec2,
                    f = NULL,
                    mel = FALSE,
                    plot = TRUE,
                    plotval = TRUE,
                    method = "spearman",
                    col = "black",
                    colval = "red",
                    cexval = 1,
                    fontval = 1,
                    xlab = NULL,
                    ylab = "Coefficient of correlation (r)",
                    type ="l",
                    ...
                    )

{
    options(warn=-1) ## to remove warning messages from cor.test()
    if(is.matrix(spec1) && ncol(spec1)==2) s1<-spec1[,2] else s1 <- spec1
    if(is.matrix(spec2) && ncol(spec2)==2) s2<-spec2[,2] else s2 <- spec2

    n1<-length(s1)
    n2<-length(s2)

    if(n1 != n2) stop("'spec1' and 'spec2' must have the same length")
    if(any(s1 < 0) | any(s2 < 0)) stop("data does not have to be in dB")
    if(sum(s1-s2)==0) stop("'spec1' and 'spec2' are the same object")
    
    mean1<-mean(s1)
    mean2<-mean(s2)
    diffx<-s1-mean1
    r1<-numeric(n1)
    r2<-numeric(n2)

    for (i in 0:(n1-1))
        {
            r1[i]<-cor(x=s1,y=c(s2[(i+1):n2],rep(0,i)),method = method)
        }

    for (i in 0:(n1-1))
        {
            r2[i+1]<-cor(x=s1,y=c(rep(0,i),s2[1:(n2-i)]),method = method)
        }

    r2<-r2[-1]
    r<-c(rev(r1),r2)

    rmax<-max(r,na.rm=TRUE)
    offset<-which.max(r)
    if(offset<=(length(r)/2)) {offsetp<-which.max(r)} else {offsetp<-which.max(r)-1}

    if(offsetp < n1)
        {
            p<-cor.test(x=s1, y=c(s2[(n1-offsetp+1):n2],rep(0,n1-offsetp)),
                        method = method)
        }
    else
        {
            p<-cor.test(x=s1, y=c(rep(0,offsetp-n1),s2[1:(n2-(offsetp-n1))]),
                        method = method)
        }
    p<-p$p.value
    

    if (!is.null(f) & mel) {f <- 2*mel(f/2)}
    if(is.null(f))
        {
            if(is.vector(spec1) & is.vector(spec2))
                {
                    if(plot) warning("'f' is missing, frequency of maximum correlation cannot be computed - No plot produced", call.=FALSE)
                    else warning("'f' is missing, frequency of maximum correlation cannot be computed", call.=FALSE)
                    corr <- list(r = r, rmax = rmax, p = p, f = NA)
                    return(corr)
                }
            else
                {
                    if(is.matrix(spec1)) f<-spec1[nrow(spec1),1]*2000*nrow(spec1)/(nrow(spec1)-1)
                    else if(is.matrix(spec2)) f<-spec2[nrow(spec2),1]*2000*nrow(spec2)/(nrow(spec2)-1)
                }
        }

    range<-c(0,f/2000)
    offsetf<-((range[2]-range[1])*(offsetp-n1))/n1
    unname(offsetf)   # to remove an x appearing as a name
    X<-seq(-range[2],range[2],length.out=2*n1-1)

    corr<-list(r = cbind(X,r), rmax=rmax, p=p, f=offsetf)

    options(warn=0) ## to reset warning messages

    if(plot)
        {
         if (mel) {w <- " kmel"} else {w <- " kHz"}
         if(is.null(xlab)) {if (mel) xlab <- "Frequency (kmel)" else xlab <- "Frequency (kHz)"}
         plot(x = X, y = r, xlab = xlab, ylab = ylab, col = col, type = type,...)
         if(plotval)
             {
                 mtext(paste(
                     "rmax = ", as.character(round(rmax,2)),
                     ", offset = ", as.character(round(offsetf,2)), w, sep=" "),
                       side=3, line=-2, col=colval, cex=cexval, font=fontval)
                 segments(x0=offsetf,y0=min(r), x1=offsetf, y1=rmax, col=colval, lty=2)
                 segments(x0=-range[2],y0=rmax, x1=offsetf, y1=rmax, col=colval, lty=2)
                 points(x=offsetf,y=rmax,pch=19,cex=1, col=colval)
             }
      return(invisible(corr))    # 2016-12-01 does not work if not wrapped by return
        }
    else  
        {
            return(corr)
        }
}


################################################################################
##                                COVSPECTRO                                         
################################################################################


covspectro<-function(
                     wave1,
                     wave2,
                     f,
                     channel = c(1,1),
                     wl = 512,
                     wn = "hanning",
                     n,
                     plot = TRUE,
                     plotval = TRUE,
                     method = "spearman",
                     col = "black",
                     colval = "red",
                     cexval = 1,
                     fontval = 1,
                     xlab = "Time (s)",
                     ylab = "Normalised covariance (cov)",
                     type ="l",
                     pb = FALSE,
                     ...
                     )

{
  input1<-inputw(wave=wave1,f=f,channel=channel[1]) ; wave1<-input1$w ; f<-input1$f ; rm(input1)
  wave2<-inputw(wave=wave2,f=f,channel=channel[2])$w

  if(n>21)
    {
      cat("please wait...")
      if(.Platform$OS.type == "windows") flush.console()
    }
  
  if(n %% 2 != 1) stop("'n' must be odd")

  n<-(n%/%2)+1

  n1<-nrow(wave1)
  n2<-nrow(wave2)

  if(n1 != n2) stop("'wave 1' and 'wave 2' must have the same length")

  step1<-seq(1,n1-wl,wl); lstep1<-length(step1)
  step2<-round(seq(1,n2,length.out=n))

                                        # wave not time shifted
  spectro1<-sspectro(wave=wave1,f=f,wl=wl,wn=wn)
  WL<-wl%/%2
  spectro2a<-array(numeric(WL*lstep1*n),dim=c(WL,lstep1,n))
  spectro2b<-array(numeric(WL*lstep1*n),dim=c(WL,lstep1,n))
  cov1<-numeric(n)
  cov2<-numeric(n)

                                        # wave time shifted
                                        # successive spectrograms
                                        # covariance of spectrogram1/spectra1 with spectrogram2/spectra1,
                                        # spectrogram1/spectra2 with spectrogram2/spectra2 and so on
                                        # diagonal of the cov matrix and mean of this diagonal
                                        # one mean cov for comparaison between 2 spectrograms for(i in step2)
  if(pb) {pbar <- txtProgressBar(min=0, max=n2, style = 3)}
  for (i in step2)
    {
      spectro2a[,,which(step2==i)]<-sspectro(wave=as.matrix(c(wave2[i:n2],rep(0,i-1))),f=f,wl=wl,wn=wn)
      spectro2a<-ifelse(spectro2a=="NaN",yes=0,no=spectro2a)
      cov1[which(step2==i)]<-mean(diag(cov(x=spectro1,y=spectro2a[,,which(step2==i)],method = method)))
      spectro2b[,,which(step2==i)]<-sspectro(wave=as.matrix(c(rep(0,i),wave2[1:(n2-i)])),f=f,wl=wl,wn=wn)
      spectro2b<-ifelse(spectro2b=="NaN",yes=0,no=spectro2b)
      cov2[which(step2==i)]<-mean(diag(cov(x=spectro1,y=spectro2b[,,which(step2==i)],method = method)))
      if(pb) {setTxtProgressBar(pbar, i)}
    }

                                        # to normalise the covariance we need covmax that is the autocovariance of spectro1
  covmax<-mean(diag(cov(x=spectro1,y=spectro1,method = method)))

                                        # discard the first value of cov2 that is already computed in cov1
  cov2<-cov2[-1]
  cov3<-c(rev(cov1),cov2)
  cov4<-cov3/covmax
  cov4max<-max(cov4)
  offset<-which.max(cov4)
  offsetp<-which.max(cov4)
  offsett<-(((offsetp*n1)/n)-n1)/f
  covar<-list(cov = cov4, covmax = cov4max, t = offsett)

  if(plot)
    {
      x<-seq(-n1/f,n1/f,length.out=(2*n)-1)
      plot(x = x, y = cov4, xlab = xlab, ylab = ylab, col = col, type = type,...)
      if(plotval)
        {
          mtext(paste(
                      "covmax = ", as.character(round(cov4max,2)),
                      ", offset = ", as.character(round(offsett,3)), "s", sep=" "),
                side=3, line=-2, col=colval, cex=cexval, font=fontval)
          segments(x0=round(offsett,3),y0=min(cov4), x1=round(offsett,3), y1=cov4max, col=colval, lty=2)
          segments(x0=-n1/f,y0=cov4max, x1=round(offsett,3), y1=cov4max, col=colval, lty=2)
          points(x=round(offsett,3),y=cov4max,pch=19,cex=1, col=colval)
        }
     return(invisible(covar))    # 2016-12-01 does not work if not wrapped by return 
    }

  else
    {
      return(covar)
    }
  if(pb) close(pbar)
}


################################################################################
##                               CREST
################################################################################

crest <- function(
                  wave,
                  f,
                  channel = 1,
                  plot = FALSE,
                  col = 2,
                  cex = 3,
                  symbol = "*",
                  ...
                  )

{  
  input<-inputw(wave=wave,f=f, channel=channel) ; wave<-input$w ; f<-input$f ; rm(input)
  max <- max(abs(wave))
  loc.max <- which(wave == max(wave))/f
  c <- max/rms(wave)

  if(plot)
    {
      if(which.max(abs(range(wave)))==1) {max2 <- -max} else {max2 <- max}# when the max is actually a negative value
      oscillo(wave=wave, f=f,...)
      text(x=loc.max, y=max2, labels=symbol, cex=cex, col=col)
    }
  return(list(C=c, val=max, loc=loc.max))
}


################################################################################
##                                CSH
################################################################################


csh<-function(
              wave,
              f,
              channel = 1,
              wl = 512,
              wn = "hanning",
              ovlp = 0,
              fftw = FALSE,
              threshold = NULL,
              plot = TRUE,
              xlab = "Times (s)",
              ylab = "Spectral Entropy",
              ylim = c(0,1.1),
              type = "l",
              ...)

{
  input<-inputw(wave=wave,f=f,channel=channel) ; wave<-input$w ; f<-input$f ; rm(input)

                                        # threshold
  if(!is.null(threshold)) wave<-afilter(wave=wave,f=f,threshold=threshold,plot=FALSE)

                                        # STFT
  n<-nrow(wave)
  step<-seq(1,n+1-wl,wl-(ovlp*wl/100)) # +1 added @ 2017-04-20
  z <- stdft(wave=wave, f=f, wl=wl, zp=0, step=step, wn=wn, fftw=fftw)
                                        # sh applied to the Fourier matrix
  h<-apply(z, MARGIN=2, FUN=sh)

  t<-seq(0,n/f,length.out=length(step))

                                        # graphic
  if(plot)
    {
      plot(x=t, y=h,
           xaxs = "i", xlab = xlab,
           yaxs = "i", ylab = ylab, ylim = ylim,
           type = type,
           ...)
      invisible(cbind(t,h))
    }
  else
    return(cbind(t,h))
}


################################################################################
##                                CUTSPEC                                       
################################################################################

cutspec <- function(spec,
                     f = NULL,
                     flim,
                     mel = FALSE,
                     norm = FALSE,
                     PMF = FALSE
                     )

{
    if(norm & PMF) stop("'norm' and 'PMF' should not be both set to TRUE")

    if(is.vector(spec))
        {
            if(is.null(f)) stop("'f' is missing and is necessary when 'spec' is a vector")
            if(mel) {
                f <- 2*mel(f/2)
		flim <- mel(flim*1000)/1000
            } 
            wl <- length(spec)*2
            if(flim[2]==f/2000) {i <- 1; j=0} else {i <- j <- 1}        ## to avoid to search over (Nyquist-1) and because first value at '0 Hz'
            flim.index.inf <- round((flim[1]*1000*wl/2)/(f/2 - f/wl)+i)
            flim.index.sup <- round((flim[2]*1000*wl/2)/(f/2 - f/wl)+j)
            specut <- spec[flim.index.inf:flim.index.sup]
            if(norm) {specut <- specut/max(specut)}
            if(PMF)  {specut <- specut/sum(specut)}     
        }

    else if(is.matrix(spec))
        {
            if(ncol(spec)>2) {stop("'spec' should not have more than two columns")}
            if(is.null(f)) {
                f <- spec[nrow(spec),1]*2000*nrow(spec)/(nrow(spec)-1)
                if (mel) {flim <- mel(flim*1000)}
            }
            wl <- nrow(spec)*2
            if(flim[2]>=(f/2-f/wl)/1000) {
                if(flim[2]>=(f/2-f/wl)/1000) {flim[2] <- (f/2-f/wl)/1000}
                i <- 1; j=0
            }
            else {i <- j <- 1}           ## to avoid to search over (Nyquist-1) and because first value at '0 Hz'
            flim.index.inf <- round((flim[1]*1000*wl/2)/(f/2 - f/wl)+i)
            flim.index.sup <- round((flim[2]*1000*wl/2)/(f/2 - f/wl)+j)
            specut <- spec[flim.index.inf:flim.index.sup, , drop=FALSE] ## JS added +1 because first value at '0 Hz'
            if(norm) {specut[,2] <- specut[,2]/max(specut[,2])}
            if(PMF)  {specut[,2] <- specut[,2]/sum(specut[,2])}
        }

    return(specut)
}


################################################################################
##                                CUTW                                         
################################################################################

cutw <- function(
                 wave,
                 f,
                 channel = 1,
                 from = NULL,
                 to = NULL,
                 choose = FALSE,
                 plot = FALSE,
                 marks = TRUE,               
                 output = "matrix",
                 ...)

{
  input <- inputw(wave=wave,f=f,channel=channel) ; wave <- input$w ; f <- input$f ; bit <- input$bit ; rm(input)

  if(choose)
    { 
      cat("choose start and end positions on the wave\n")
      if(.Platform$OS.type == "windows") flush.console()
      oscillo(wave,f=f)
      coord <- locator(n=2)
      from <- coord$x[1] ; a <- round(from*f)+1 ; abline(v=from,col=2,lty=2)
      to <- coord$x[2]   ; b <- round(to*f)     ; abline(v=to,col=2,lty=2)
    }
  else if(!is.null(from)|!is.null(to))
    {
      if(is.null(from) && !is.null(to)) {a <- 1 ; b <- round(to*f)}
      if(!is.null(from) && is.null(to)) {a <- round(from*f)+1 ; b <- length(wave)}
      if(!is.null(from) && !is.null(to))
        {
          if(from>to) stop("'from' cannot be superior to 'to'")
          if(from==0) {a <- 1} else a <- round(from*f)+1
          b <- round(to*f)
        }
    }

  wavecut <- as.matrix(wave[a:b,])

  wavecut <- outputw(wave=wavecut, f=f, bit=bit, format=output)

  if(plot)
    {
      def.par <- par(no.readonly = TRUE)
      on.exit(par(def.par))
      par(mfrow=c(2,1),oma=c(0,0.1,0,0))
      oscillo(wave,f=f,...)
      title(main="original")
      if(marks)
        {
          abline(v=from, col="red", lty=2)
          abline(v=to, col="red", lty=2)
        }
      oscillo(wavecut,f=f,...)
      title(main="cut")
      invisible(wavecut)
    }
  else return(wavecut)
}


################################################################################
##                                DBSCALE                                        
################################################################################

dBscale<-function
(
 collevels,
 palette = spectro.colors,
 side = 4,
 textlab = "Amplitude\n(dB)",
 cexlab = 0.75,
 fontlab = 1,
 collab = "black",
 colaxis = "black",
 ...
 )

{
  plot.new()
  levels<-collevels
  col <- palette(length(collevels) - 1)
  par(las=1)
  
  if(side == 2 | side == 4)
    {    
      plot.window(xlim = c(0, 1), ylim = range(collevels), xaxs = "i",
                  yaxs = "i")
      mtext(textlab, side=3, outer=FALSE, line=1.5, adj=0, font=fontlab, cex=cexlab, col=collab)
      rect(xleft=0, ybottom=levels[-length(levels)], xright=0.95, ytop=levels[-1],
           col = col, lty=0, border = TRUE)
      segments(x0=0,y0=max(collevels),x1=0.95,y1=max(collevels),col=colaxis)
      segments(x0=0,y0=min(collevels),x1=0.95,y1=min(collevels),col=colaxis)          
      abline(v=c(0,0.95),col=colaxis)
      if(side == 2) axis(2,col=colaxis,col.axis=colaxis,...)
      if(side == 4) axis(4,pos=0.95,col=colaxis,col.axis=colaxis,...)
    }

  if(side == 1  | side == 3)
    {    
      plot.window(xlim = range(collevels), ylim = c(0, 1), xaxs = "i",
                  yaxs = "i")
      mtext(textlab, side=3, outer=FALSE, line=1.5, adj=0, font=fontlab, cex=cexlab, col=collab)
      rect(xleft=levels[-length(levels)], ybottom=0, xright=levels[-1], ytop=0.95, col = col, lty=0)
      segments(x0=min(collevels),y0=0,x1=min(collevels),y1=0.95,col=colaxis)
      segments(x0=max(collevels),y0=0,x1=max(collevels),y1=0.95,col=colaxis)       
      abline(h=c(0,0.95),col=colaxis)
      if(side == 1) axis(1,col=colaxis,col.axis=colaxis,...)
      if(side == 3) axis(3,pos=0.95,col=colaxis,col.axis=colaxis,...)
    }    
}    


################################################################################
##                                DBWEIGHT
################################################################################

dBweight <- function(
                    f,
                    dBref=NULL
                    )
  
{
  # dBA
  num <- (12200^2*f^4)
  den <- (f^2 + 20.6^2)* sqrt((f^2 + 107.7^2)*(f^2 + 737.9^2))*(f^2 + 12200^2)
  A <- 2 + 20*log10(num/den)
  A <- ifelse(is.infinite(A), yes=NA, no=A)

  # dBB
  num <- (12200^2*f^3);
  den <- (f^2 + 20.6^2)*sqrt((f^2 + 158.5^2))*(f^2 + 12200^2)
  B <- 0.17 + 20*log10(num/den)
  B <- ifelse(is.infinite(B), yes=NA, no=B)

  # dBC
  num <- (12200^2*f^2);
  den <- (f^2 + 20.6^2)*(f^2 + 12200^2)
  C <- 0.06 + 20*log10(num/den)
  C <- ifelse(is.infinite(C), yes=NA, no=C)

  # dBD
  a <- f/(6.8966888496476*10^-5);
  h <- ((1037918.48 - f^2)^2 + (1080768.16*f^2))/((9837328 - f^2)^2 + (11723776*f^2))
  b <- sqrt(h/((f^2 + 79919.29)*(f^2 + 1345600)))
  D <- 20*log10(a*b)
  D <- ifelse(is.infinite(D), yes=NA, no=D)

  # dB ITU-R 468
  # added @ 20108-06-18 by Andrey Anikin
  # https://en.wikipedia.org/wiki/ITU-R_468_noise_weighting
  h1 <- -4.737338981378384 * 10^(-24) * f^6 + 
      2.04382833606125 * 10^(-15) * f^4 - 
      1.363894795463638 * 10^(-7) * f^2 + 1
  h2 <- 1.306612257412824 * 10^(-19) * f^5 - 
      2.118150887518656 * 10^(-11) * f^3 + 
      5.559488023498642 * 10^(-4) * f
  R_ITU <- 1.246332637532143 * 10^(-4) * f / sqrt(h1^2 + h2^2)
  ITU <- 18.2 + 20 * log10(R_ITU)

  # result
  result <- list(A=A, B=B, C=C, D=D, ITU=ITU)
  if(!is.null(dBref)) result <- list(A=dBref+A, B=dBref+B, C=dBref+C, D=dBref+D, ITU=dBref+ITU)
  return(result)
}


################################################################################
##                                DELETEW
################################################################################

deletew<-function(
                  wave,
                  f,
                  channel = 1,
                  from = NULL,
                  to = NULL,
                  choose = FALSE,
                  plot = FALSE,
                  marks = TRUE,
                  output = "matrix",
                  ...)

{
  input<-inputw(wave=wave,f=f,channel=channel) ; wave<-input$w ; f<-input$f ; bit <- input$bit ; rm(input)

  if(choose)
    { 
      cat("choose start and end positions on the wave\n")
      if(.Platform$OS.type == "windows") flush.console()
      oscillo(wave,f=f)
      coord<-locator(n=2)
      from<-coord$x[1]; a<-round(from*f)+1 ; abline(v=from,col=2,lty=2)
      to<-coord$x[2];   b<-round(to*f)     ; abline(v=to,col=2,lty=2)
    }
  else if(!is.null(from)|!is.null(to))
    {
      if(is.null(from) && !is.null(to)) {a<-1; b<-round(to*f)}
      if(!is.null(from) && is.null(to)) {a<-round(from*f)+1 ; b<-length(wave)}
      if(!is.null(from) && !is.null(to))
        {
          if(from>to) stop("'from' cannot be superior to 'to'")
          if(from==0) {a<-1} else a<-round(from*f)+1
          b<-round(to*f)
        }
    }
  wavecut <- wave[-(a:b),]

  wavecut <- outputw(wave=wavecut, f=f, format=output)

  if(plot)
    {
      def.par <- par(no.readonly = TRUE)
      on.exit(par(def.par))
      par(mfrow=c(2,1))
      oscillo(wave,f=f,k=1,j=1,...)
      title(main="original")
      if(marks)
        {
          abline(v=from, col="red", lty=2)
          abline(v=to, col="red", lty=2)
        }
      oscillo(wavecut,f=f,k=1,j=1,...)
      title(main="after deletion")
      invisible(wavecut)
    }
  else return(wavecut)
}


################################################################################
## DFREQ                                         
################################################################################

dfreq <- function(
                wave,
                f,
                channel = 1,
                wl = 512,
                wn = "hanning",
                ovlp = 0,
                fftw = FALSE,
                at = NULL,
                tlim = NULL,
                threshold = NULL,
                bandpass = NULL,
                clip = NULL,  
                plot = TRUE,
                xlab = "Times (s)",
                ylab = "Frequency (kHz)",
                ylim = c(0,f/2000),
                ...)

{
  # error messages
  if(!is.null(at) && ovlp != 0) stop("The 'ovlp' argument cannot bue used in conjunction with the arguement 'at'.")
  if(!is.null(clip)) {if(clip <=0 | clip >= 1) stop("'clip' value has to be superior to 0 and inferior to 1")}

  # input
  input<-inputw(wave=wave,f=f,channel=channel) ; wave<-input$w ; f<-input$f ; rm(input)

  # time limits
  if(!is.null(tlim)) wave<-cutw(wave,f=f,from=tlim[1],to=tlim[2])                                      
  
  # amplitude threshold
  if(!is.null(threshold)) {wave<-afilter(wave=wave,f=f,threshold=threshold,plot=FALSE)}

  # Position(s)
  n<-nrow(wave)
  if(!is.null(at))
       {
       step <- at*f
       N <- length(step)
       if(step[1]<=0) {step[1] <- 1}
       if(step[N]+(wl/2) >= n) {step[N] <- n-wl}
       x <- c(0, at, n/f)
       }
  else {
       step <- seq(1,n+1-wl,wl-(ovlp*wl/100)) # +1 added @ 2017-04-20
       N <- length(step)
       x <- seq(0, n/f, length.out=N)
       }

  # Fourier
  step <- round(step)
  y1 <- stdft(wave=wave, f=f, wl=wl, zp=0, step=step, wn=wn, fftw=fftw) 
  
  # bandpass filter, values outside the bandpass limits are replaced by 0 values
  if(!is.null(bandpass))
    {
      if(length(bandpass)!=2) stop("'The argument 'bandpass' should be a numeric vector of length 2'")
      if(bandpass[1] > bandpass[2]) stop("The first element of 'bandpass' has to be inferior to the second element, i.e. bandpass[1] < bandpass[2]")
      if(bandpass[1] == bandpass[2]) stop("The limits of the bandpass have to be different")
      lowlimit <-round((wl*bandpass[1])/f)
      upperlimit <-round((wl*bandpass[2])/f)
      y1[-(lowlimit:upperlimit),] <- 0
   }

  # Maximum search
  maxi <- apply(y1, MARGIN=2, FUN=max)
  y2 <- apply(y1, MARGIN=2, FUN=which.max)
  y2[which(maxi==0)] <- NA
 
  # discards peaks with an amplitude lower than the clip value
  if(!is.null(clip))
    {
      maxi <- apply(y1, MARGIN=2, FUN=max)
      y2[which(maxi < clip)] <- NA
    }
  # converts into frequency
  y <- (f*y2)/(1000*wl) - f/(1000*wl)

  if(!is.null(at)) {y <- c(NA, y, NA)}

  if(plot)
    {
      plot(x=x, y=y,
           xaxs="i", xlab = xlab,
           yaxs="i", ylab = ylab, ylim = ylim,
           ...)
      invisible(cbind(x,y))      
    }
  else
    return(cbind(x,y))
}


################################################################################
##                                DIFFCUMSPEC                                        
################################################################################

diffcumspec <- function(spec1,
                        spec2,
                        f = NULL,
                        mel = FALSE,
                        plot = FALSE,
                        type = "l",
                        lty = c(1,2),
                        col = c(2,4,8),
                        flab = NULL,
                        alab = "Cumulated amplitude",
                        flim = NULL, 
                        alim = NULL,
                        title = TRUE,
                        legend = TRUE,
                        ...
                        )

    {
        leg <- c(as.character(deparse(substitute(spec1))),as.character(deparse(substitute(spec2))))

        if(is.matrix(spec1) && ncol(spec1)==2) s1 <- spec1[,2] else s1 <- spec1
        if(is.matrix(spec2) && ncol(spec2)==2) s2 <- spec2[,2] else s2 <- spec2

        n1 <- length(s1)
        n2 <- length(s2)

        if(n1 != n2) stop("spec1 and spec2 must have the same length")
        
        if(any(is.na(s1)) | any(is.na(s2)))
            {
                D <- NA
                warning("The data set contains 'NA' values. The returned values have been set to NA.", call.=FALSE)
                return(D)
            } 
        else
            {
                if(any(s1 < 0) | any(s2 < 0)) stop("spectra (spec 1 and/or spec 2) do not have to be in dB")
                if(sum(s1)==0) {warning(paste("Caution!, spec1 is a null spectrum"), call.=FALSE)}
                if(sum(s2)==0) {warning(paste("Caution!, spec2 is a null spectrum"), call.=FALSE)}

		s1<-s1/sum(s1)
  		s2<-s2/sum(s2)
                cum.s1 <- cumsum(s1)
                cum.s2 <- cumsum(s2)     
                D <- sum(abs(cum.s1 - cum.s2))/n1

                if(plot)
                    {
                        if (!is.null(f) & mel) {f <- 2*mel(f/2)}
                        if(is.null(f))
                            {
                                if(is.vector(spec1) & is.vector(spec2)) stop("'f' is missing")  
                                else
                                    {
                                        if(is.matrix(spec1)) f<-spec1[nrow(spec1),1]*2000*nrow(spec1)/(nrow(spec1)-1)
                                        else if(is.matrix(spec2)) f<-spec2[nrow(spec2),1]*2000*nrow(spec2)/(nrow(spec2)-1)
                                    }
                            }
			x  <- seq(0, f/2*(n1-1)/n1, length.out=n1)/1000
                        st <- 0
                        en <- (f/2*(n1-1)/n1)/1000
      			if(is.null(alim)) alim<-c(0,1)
      			if(is.null(flim)) flim<-c(0,f/2000)
			if(is.null(flab)) {if (mel) flab <- "Frequency (kmel)" else flab <- "Frequency (kHz)"}
                        plot(x=x, y=cum.s1,
                             col=col[1], lty=lty[1], type = "n",
                             xlim = flim, xaxs = "i", xlab = flab, 
                             ylim = alim, yaxs = "i", ylab = alab, ...)
                        polygon(x=c(seq(st,en,length.out=n1),seq(en,st,length.out=n1)),
                                y=c(cum.s1,rev(cum.s2)),
                                col=col[3], border=NA)
                        lines(x=x, y=cum.s1, type=type, lty=lty[1], col=col[1])
                        lines(x=x, y=cum.s2, type=type, lty=lty[2], col=col[2])
                        if(legend) legend("topleft", col=c(col[1],col[2]) , lty=c(lty[1], lty[2]), legend=leg, bty="n")
                        if(title) title(main=paste("D = ", round(D, 3)))
                        invisible(D)
                    }
                else return(D)
            }
    }


################################################################################
##                                DIFFENV                                        
################################################################################

diffenv <- function(
                  wave1,
                  wave2,
                  f,
                  channel = c(1,1),
                  envt = "hil",
                  msmooth = NULL,
                  ksmooth = NULL,
                  plot = FALSE,
                  lty1 = 1,
                  lty2 = 2,
                  col1 = 2,
                  col2 = 4,
                  cold = 8,
                  xlab = "Time (s)",
                  ylab = "Amplitude",
                  ylim = NULL,
                  legend = TRUE,
                  ...
                  )

{
  if(length(channel)!=2) stop("'channel' should be a numeric vector of length 2")
    
  leg<-c(as.character(deparse(substitute(wave1))),as.character(deparse(substitute(wave2))))

  input1<-inputw(wave=wave1,f=f,channel=channel[1]) ; wave1<-input1$w ; f<-input1$f ; rm(input1)
  wave2<-inputw(wave=wave2,f=f,channel=channel[2])$w

  env1<-env(wave=wave1,f=f,envt=envt,msmooth=msmooth,ksmooth=ksmooth,plot=FALSE)
  env2<-env(wave=wave2,f=f,envt=envt,msmooth=msmooth,ksmooth=ksmooth,plot=FALSE)

  n1<-length(env1)
  n2<-length(env2)
  if(n1 != n2) stop("wave1 and wave2 should have the same length")

  if(!is.null(msmooth)|!is.null(ksmooth)) {f<-f*n1/nrow(wave1)}

  envPMF1<-env1/sum(env1)
  envPMF2<-env2/sum(env2)

  denv<-sum(abs(envPMF1-envPMF2))/2

  if(plot)
    {
      x<-seq(0,n1/f,length.out=n1)
      if(is.null(ylim)) ylim<-c(0,max(envPMF1,envPMF2))
      plot(x=x, y=envPMF1, type="n",
           xaxs="i", xlab = xlab, 
           ylim=ylim, yaxs="i", ylab=ylab, ...)
      polygon(x=c(x,rev(x)),
              y=c(envPMF1,rev(envPMF2)),
              col=cold,
              border=NA)
      lines(x, envPMF1, lty=lty1, col=col1)
      lines(x, envPMF2, lty=lty2, col=col2)
      box()
      if(legend) {legend("topleft", col=c(col1,col2),lty=c(lty1,lty2),legend=leg)}
      invisible(denv)
    }
  return(denv)
}


################################################################################
##                                DIFFSPEC                                         
################################################################################

diffspec <- function(spec1,
                     spec2,
                     f = NULL,
                     mel = FALSE,
                     plot = FALSE,
                     type = "l",
                     lty = c(1,2),
                     col = c(2,4,8),
                     flab = NULL,
                     alab = "Amplitude",
                     flim = NULL,
                     alim = NULL,
                     title = TRUE,
                     legend = TRUE,
                     ...
                     )

{
    leg<-c(as.character(deparse(substitute(spec1))),as.character(deparse(substitute(spec2))))

    if(is.matrix(spec1) && ncol(spec1)==2) s1<-spec1[,2] else s1 <- spec1
    if(is.matrix(spec2) && ncol(spec2)==2) s2<-spec2[,2] else s2 <- spec2

    n1<-length(s1)
    n2<-length(s2)

    if(n1 != n2) stop("spec1 and spec2 must have the same length")

    if(any(is.na(s1)) | any(is.na(s2)))
        {
            dspec <- NA
            warning("The data set contains 'NA' values. The returned values have been set to NA.", call.=FALSE)
            return(dspec)
        } 
    else
	{
            if(any(s1 < 0) | any(s2 < 0)) stop("spectra (spec 1 and/or spec 2) do not have to be in dB")
            if(sum(s1)==0) {warning(paste("Caution!, spec1 is a null spectrum"), call.=FALSE)}
            if(sum(s2)==0) {warning(paste("Caution!, spec2 is a null spectrum"), call.=FALSE)}

            s1<-s1/sum(s1)
            s2<-s2/sum(s2)

            dspec<-sum(abs(s1-s2))/2

            if(plot)
                {
                    if (!is.null(f) & mel) {f <- 2*mel(f/2)}
                    if(is.null(f))
                        {
                            if(is.vector(spec1) & is.vector(spec2)) stop("'f' is missing")  
                            else
                                {
                                    if(is.matrix(spec1)) f<-spec1[nrow(spec1),1]*2000*nrow(spec1)/(nrow(spec1)-1)
                                    else if(is.matrix(spec2)) f<-spec2[nrow(spec2),1]*2000*nrow(spec2)/(nrow(spec2)-1)
                                }
                        }
                    x<-seq(0, f/2*(n1-1)/n1, length.out=n1)/1000
                    st <- 0
                    en <- (f/2*(n1-1)/n1)/1000
                    if(is.null(alim)) alim<-c(0,max(s1,s2))
                    if(is.null(flim)) flim<-c(0,f/2000)
                    if(is.null(flab)) {if (mel) flab <- "Frequency (kmel)" else flab <- "Frequency (kHz)"}
                    plot(x=x, y=s1, type="n",
                         xlim=flim, xlab=flab, xaxs="i",
                         ylim=alim, ylab=alab, yaxs="i",...)
                    polygon(x=c(seq(st,en,length.out=n1),seq(en,st,length.out=n1)),
                            y=c(s1,rev(s2)),
                            col=col[3],
                            border=NA)
                    lines(x=x, y=s1, type=type, lty=lty[1], col=col[1])
                    lines(x=x, y=s2, type=type, lty=lty[2], col=col[2])
                    box()
                    if(legend) legend("topleft", col=c(col[1],col[2]),lty=c(lty[1],lty[2]),legend=leg, bty="n")
                    if(title) title(main=paste("D = ", round(dspec, 3)))
                    invisible(dspec)
                }
            else return(dspec)
        }
}



################################################################################
##                               DIFFWAVE                                        
################################################################################

diffwave<-function(
                   wave1,
                   wave2,
                   f,
                   channel = c(1,1),
                   wl = 512,
                   envt= "hil",
                   msmooth = NULL,
                   ksmooth = NULL
                   )

{
  input1<-inputw(wave=wave1,f=f,channel=channel[1]) ; wave1<-input1$w ; f<-input1$f ; rm(input1)
  wave2<-inputw(wave=wave2,f=f,channel=channel[2])$w

                                        # spectral difference
  spec1<-meanspec(wave=wave1,f=f,wl=wl,PMF=TRUE,plot=FALSE)
  spec2<-meanspec(wave=wave2,f=f,wl=wl,PMF=TRUE,plot=FALSE)
  DF<-diffspec(spec1=spec1,spec2=spec2,f=f,plot=FALSE)

                                        # temporal difference 
  DE<-diffenv(wave1=wave1,wave2=wave2,f=f,msmooth=msmooth,ksmooth=ksmooth,plot=FALSE)

  z<-DF*DE
  return(z)
}


################################################################################
##                               DISCRETS                                        
################################################################################

discrets <- function (x,
                      symb = 5,
                      collapse = TRUE,
                      plateau = 1
                      ) 

{
    ## DATA INPUT
    if (symb != 3 && symb != 5) 
        stop("'symb' should be set to 3 or 5")
    x <- inputw(wave = x, f = NULL)$w
    n <- length(x)
    ## 5 SYMBOLS
    if(symb==5)
        {
            ## from the second point to the n-1 point
            s <- character(n-2)
            if(plateau == 1){     ## Cazelles et al definition, a plateau is encoded as PFP      
                for (i in 1:(n-2)){
                    if(x[i]<=x[i+1]   & x[i+1]<x[i+2])  s[i+1]<-"I"  # increase
                    if(x[i]<=x[i+2]   & x[i+2]<=x[i+1]) s[i+1]<-"P"  # peak
                    if(x[i+1]<x[i]    & x[i]<=x[i+2])   s[i+1]<-"T"  # trough
                    if(x[i+1]<x[i+2]  & x[i+2]<=x[i])   s[i+1]<-"T"  # trough
                    if(x[i+2]<x[i]    & x[i]<=x[i+1])   s[i+1]<-"P"  # peak
                    if(x[i+2]<=x[i+1] & x[i+1]<x[i])    s[i+1]<-"D"  # decrease
                    if(x[i]==x[i+1]   & x[i+1]==x[i+2]) s[i+1]<-"F"  # flat
                }
            }
            else if (plateau == 2) { ## Lellouch definition, a plateau is encoded as IFD
                for (i in 1:(n - 2)) {
                    if (x[i] < x[i + 1] & x[i + 1] > x[i + 2])   s[i + 1] <- "P" # peak
                    if (x[i] > x[i + 1] & x[i + 1] < x[i + 2])   s[i + 1] <- "T" # trough
                    if (x[i] <= x[i + 1] & x[i + 1] <= x[i + 2]) s[i + 1] <- "I" # increase
                    if (x[i] >= x[i + 1] & x[i + 1] >= x[i + 2]) s[i + 1] <- "D" # decrease
                    if (x[i] == x[i + 1] & x[i + 1] == x[i + 2]) s[i + 1] <- "F" # flat
                }
            }
        }
    else if (symb == 3) {
        s <- character(n - 1)
        for (i in 1:(n - 1)) {
            if (x[i] == x[i + 1]) 
                s[i + 1] <- "F"
            if (x[i] < x[i + 1]) 
                s[i + 1] <- "I"
            if (x[i] > x[i + 1]) 
                s[i + 1] <- "D"
        }
    }
    s <- s[-1]
    if (collapse) {s <- paste(s, collapse = "")}
    return(s)
}


################################################################################
##                                DRAWENV
################################################################################

drawenv<-function(
                  wave,
                  f,
                  channel = 1, 
                  n=20,
                  plot=FALSE,
                  listen = FALSE,                  
                  output = "matrix"
                  )

{
                                        # input
  input<-inputw(wave=wave,f=f,channel=channel) ; wave<-input$w ; f<-input$f ; bit <- input$bit ; rm(input)

  wave<-wave/max(abs(wave))
  wave<-rmoffset(wave, f=f)

                                        # interactive graph
  oscillo(wave=wave,f=f)
  cat("choose points on the positive amplitude side of the wave\nto change the amplitude profile (amplitude envelope)\n")
  if(.Platform$OS.type == "windows") flush.console()
  coord<-locator(n=n,type="p",col=2)

                                        # coordinates ; ordered following x if positions are not localised in order along the x-time axis
  X<-coord$x ; x<-round(X[order(X)]*f)
  Y<-coord$y ; y<-Y[order(X)]
  if(any(X<0)) stop("point localization cannot be on the negative part of the time axis")
  if(any(X>(nrow(wave)/f))) stop("point localization cannot be outside the positive part of the time axis")
  if(any(Y<0)) stop("point localization cannot be on the negative part of the amplitude axis")

                                        # profile generation
  profile<-numeric(nrow(wave))
  profile[1:x[1]]<-seq(0,y[1],length.out=x[1])
  for(i in 1:(length(x)-1))
    {
      profile[x[i]:x[i+1]]<-seq(y[i],y[i+1],length.out=x[i+1]-x[i]+1)
    }
  profile[x[length(x)]:length(profile)]<-seq(y[length(x)],0,length.out=length(profile)-x[length(x)]+1)

                                        # new wave generation
  wave2<-rmam(wave,f=f)
  wave2<-wave2/max(abs(wave2))
  wave3<-wave2[,1]*profile

  wave3 <- outputw(wave=wave3, f=f, format=output)

  if(plot)
    {
      dev.new()
      oscillo(wave3,f=f)
      if(listen) {listen(wave3,f=f)}
      invisible(wave3)
    }
  else
    {
      if(listen) {listen(wave3,f=f)}
      return(wave3)
    }
}


################################################################################
##                               DRAWFILTER
################################################################################

drawfilter <- function(f,   
                       n = 256,
                       continuous = TRUE,
                       discrete = TRUE
                       )
    
{
    if (!(continuous || discrete)) {stop("continuous or discrete must be TRUE")}
    f <- f/2000 
    plot.new()
    par(xaxs="i", yaxs="i", las=1)  
    plot(x=NA, xlim=c(0,f), ylim=c(0,1), type="n", xlab="Frequency (kHz)", ylab="Amplitude") 
    grid() 
    axis(1)
    axis(2)
    d <- seq(0,f,length.out=n+1)[1:n]
    g <- array(0, n)
    if (continuous)
        {cat("Please draw a continuous shape\n") 
         if (.Platform$OS.type == "windows") 
             flush.console()
         a <- locator(n=n, type = "l")
         if (length(a$x) < 2) {warning ("Please choose at least two points for the continuous filter\n")}
         else
             {if (any (diff(a$x) <= 0)) {stop("Values should increase along the x (frequency) axis")} 
              e <- array(0, n)
              for (i in 1:n)
                  {if (length(which(d[i] >= a$x)) == 0)
                       {e[i] <- 0}
                   else e[i] <- max(which(d[i] > a$x))}
              for (i in 1:n)
                  {if (e[i] == 0 || e[i] == length(a$x)) {g[i] <- 0}
                   else {h <- a$y[e[i]] + (a$y[e[i]+1] - a$y[e[i]])*(d[i] - a$x[e[i]])/(a$x[e[i]+1] - a$x[e[i]])
                         g[i] <- (h+abs(h))/2}
               }
          }
     }
    if (discrete)
        {cat("Please choose discrete frequencies\n")
         if (.Platform$OS.type == "windows") 
             flush.console()
         b <- locator(n=n, type = "p")
         if (length(b$x) < 1) {warning ("Locate at least one point for the discrete filter\n")}
         j <- round(b$x*n/f)
         j[j>=-5 & j < 0] <- 0
         j[j==n] <- n-1
         k <- b$y
         p <- j[j>=0 & j<=n-1 & k>=0 & k<=1.03]
         q <- k[j>=0 & j<=n-1 & k>=0 & k<=1.03]
         ext <- length(b$x) - length(p)
         red <- 0
         m <- array(0,n)
         if (length(p) >=1)
             {for(i in 1:length(p))
                  {if (m[p[i]+1] != 0) {red <- red + 1}
                   m[p[i]+1] <- m[p[i]+1] + q[i]
               }}
         if (ext > 0) {warning (cat(as.character(ext), "Discrete points out of the frame were not taken in account\n"))}
         if (red > 0) {warning (cat(as.character(red), "Duplications in discrete frequencies were found. For such frequencies, the sum of input amplitudes was used\n"))}
         g[m>0] <- m[m>0]}
    return(cbind(freq=d,amp=g))  ## JS changed the name of the columns to make it more clear
}


################################################################################
##                               DURATION
################################################################################

duration <- function(wave, f, channel = 1){
    input <- inputw(wave=wave,f=f, channel=channel)
    return(length(input$w) / input$f)
}

################################################################################
##                               DYNOSCILLO
################################################################################

dynoscillo <- function(
                  wave,
                  f,
                  channel = 1,
                  wd = NULL,
                  wl = NULL,
                  wnb = NULL,
                  title = TRUE,
                  ...)                  

{
# STOP MESSAGES
  if(is.null(wl) && is.null(wnb) && is.null(wd)) stop("Either 'wd', 'wl' or 'wnb' has to be set.")
  if(!is.null(wl) && !is.null(wnb) && !is.null(wd)) stop("'wd', 'wl', 'wnb' cannot be used in the same time. Please choose one.")
  if(!is.null(wl) && !is.null(wnb)) stop("'wl' and 'wnb' cannot be used in the same time. Please choose one.")
  if(!is.null(wl) && !is.null(wd)) stop("'wl' and 'wd' cannot be used in the same time. Please choose one.")
  if(!is.null(wnb) && !is.null(wd)) stop("'wnb' and 'wd' cannot be used in the same time. Please choose one.")
  
# INPUT
  input<-inputw(wave=wave,f=f, channel=channel) ; wave<-input$w ; f<-input$f ; rm(input)

# SECTIONS
  n <- nrow(wave)
  if(!is.null(wd)) {wl <- wd*f}
  if(!is.null(wnb)) {wl <- round(n/wnb); wd <- wl/f}
  if(!is.null(wf)) {wd <- wl/f}
  step <- seq(from = 0, to = n-wl, by = wl)/f
  lstep <- length(step)
  pos <- 1:lstep
  
# PLOT  
      plot.osc<-function(panel)
        {
          with(panel,
                  {
                     max <- max(abs(wave))
                     soscillo(wave = wave, f = f, from = step[pos], to = step[pos]+wd, ylim=c(-max, max), tickup=max,...)
                     if(title) title(main=paste(pos,"/",lstep, "\n [wd = ", round(wd,3), " s, ", "wl = ", wl, "]", sep=""))
                  }
               )
          panel
        }

      osc.panel <- rpanel::rp.control("Position")
      rpanel::rp.slider(osc.panel, pos, from = 1, to = lstep, resolution = 1,
                title = "Position along the signal", action = plot.osc)
}


################################################################################
##                               DYNSPEC
################################################################################

dynspec<-function(
                  wave,
                  f,
                  channel = 1,
                  wl = 512,
                  wn = "hanning",
                  zp = 0,
                  ovlp = 0,
                  fftw = FALSE,
                  norm = FALSE,
                  dB = NULL,
                  dBref = NULL,
                  plot = TRUE,
                  title = TRUE,
                  osc = FALSE,
                  tlab = "Time (s)",
                  flab = "Frequency (kHz)",
                  alab = "Amplitude",
                  alim = NULL,
                  flim = c(0,f/2000),
                  type ="l",
                  from = NULL,
                  to = NULL,
                  envt = NULL,
                  msmooth = NULL,
                  ksmooth = NULL,
                  colspec = "black",
                  coltitle = "black",
                  colbg = "white",
                  colline = "black",
                  colaxis = "black",
                  collab = "black",
                  cexlab = 1,
                  fontlab = 1,
                  colwave = "black",
                  coly0 = "lightgrey",
                  colcursor = "red",
                  bty = "l"
                  )

{
## STOP MESSAGES
  if(is.logical(dB)) stop("'dB' is no more a logical. Please see the documentation: help(spec).")
  if(!is.null(dB) && all(dB!=c("max0","A","B","C","D")))
    stop("'dB' has to be one of the following character strings: 'max0', 'A', 'B', 'C' or 'D'")
  if (!is.null(envt) && osc) stop("'envt' and 'osc' cannot be used together")  ## ++
  

## INPUT
  input<-inputw(wave=wave,f=f,channel=channel) ; wave<-input$w ; f<-input$f ; rm(input)

  if(!is.null(from)|!is.null(to))
    {
      if(is.null(from) && !is.null(to)) {a<-1; b<-round(to*f)}
      if(!is.null(from) && is.null(to)) {a<-round(from*f)+1; b<-length(wave)}
      if(!is.null(from) && !is.null(to))
        {
          if(from > to) stop("'from' cannot be superior to 'to'")
          if(from==0) {a<-1} else {a<-round(from*f)+1}
          b<-round(to*f)
        }
      wave<-as.matrix(wave[a:b,])
    }

## STFT
  n <- nrow(wave)
  step <- seq(1,n+1-wl,wl-(ovlp*wl/100)) # +1 added @ 2017-04-20
  lstep <- length(step)
  z <- stdft(wave=wave, f=f, wl=wl, zp=0, step=step, wn=wn, norm=norm, fftw=fftw)
 
## FREQUENCY DATA 
  x <- seq(0, (f/2)-(f/(wl+zp)), by=f/(wl+zp))/1000

## DB
  if(!is.null(dB))
    {
      if(is.null(dBref)) z<-20*log10(z) else z<-20*log10(z/dBref)
      if(dB == "max0") z <- z
      if(dB == "A") z <- dBweight(x*1000, dBref = z)$A 
      if(dB == "B") z <- dBweight(x*1000, dBref = z)$B 
      if(dB == "C") z <- dBweight(x*1000, dBref = z)$C 
      if(dB == "D") z <- dBweight(x*1000, dBref = z)$D
    }
  
  if(plot)
    {
      if(is.null(alim))
        {
          alim<-c(0,max(z)+0.05)
          if(norm && is.null(dB)) alim<-c(0,1.1)
          if(!is.null(dB)) alim<-c(min(z),10)
        }     
      pos<-1:lstep
      poslabel<-numeric(lstep)
      for (i in 1:lstep){poslabel[i]<-round((step[i]+((step[2]-step[1])/2))/f,3)}

      plot.spec<-function(panel)
        {
          with(panel,
               {
                 par(bg=colbg, col=colline)
                 if(osc | !is.null(envt)) {layout(c(1,2),heights=c(2.5,1)); par(mar=c(4.5,4,3,2))}
                 plot(x=x,y=z[,pos],
                      xaxs = "i", xlab = flab, xlim = flim,
                      yaxs = "i", yaxt = "s", ylab = alab, ylim = alim,
                      col = colspec, col.axis = colaxis,
                      col.lab = collab, cex.lab = cexlab, font.lab = fontlab,
                      type = type, las = 1)

                 if(title)
                   {
                     nc <- ncol(z)
                     title(main=paste(pos, "/", nc, " - Position along the signal = ", poslabel[pos], "s", sep=""),
                           col.main=coltitle)
                   }
                 if(title==FALSE) title(main="")
                 if(is.character(title)) title(main = paste(title), col.main = coltitle)

                 if(osc)
                   {
                     par(mar=c(4.5,4,0.5,2))
                     max <- max(abs(wave))
                     soscillo(wave = wave, f = f,
                              colwave = colwave, collab = collab,
                              cexlab = cexlab, fontlab = fontlab, colline = colline,
                              colaxis = colaxis, coly0 = coly0, tlab=tlab, bty = bty,
                              tickup=max(abs(wave),na.rm=TRUE),
                              ylim=c(-max,max))
                     abline(v=poslabel[pos], col=colcursor)
                   }
                 else if(!is.null(envt))
                   {
                     par(mar=c(4.5,4,0.5,2))
                     env(wave = wave, f = f, k=1, j=1,
                         envt = envt, msmooth = msmooth, ksmooth = ksmooth,
                         colwave = colwave, collab = collab,
                         cexlab = cexlab, fontlab = fontlab, colline = colline,
                         colaxis = colaxis, coly0 = coly0, bty = bty)
                     abline(v=poslabel[pos], col=colcursor)
                   }
               }
               )
          panel
        }

      spec.panel <- rpanel::rp.control("Position")
      rpanel::rp.slider(spec.panel, pos, from=1, to=lstep, resolution=1,
                title = "Position along the signal", action=plot.spec)
    }
  else return(list(time=seq(0,n/f,length.out=length(step)), freq=x, amp=z))
}


################################################################################
##                               DYNSPECTRO
################################################################################

dynspectro <- function(
                       wave,
                       f,
                       channel = 1,
                       slidframe = 10,
                       wl = 512,
                       wn = "hanning",
                       zp = 0,
                       ovlp = 75,
                       fftw = FALSE,
                       dB = TRUE, 
                       plot = TRUE,
                       title = TRUE,
                       osc = FALSE,
                       tlab = "Time (s)",
                       flab = "Frequency (kHz)", 
                       alab = "Amplitude",
                       from = NULL,
                       to = NULL,
                       collevels = NULL, 
                       palette = spectro.colors,
                       envt = NULL,
                       msmooth = NULL,
                       ksmooth = NULL, 
                       coltitle = "black",
                       colbg = "white",
                       colline = "black", 
                       colaxis = "black",  
                       collab = "black", 
                       cexlab = 1,  
                       fontlab = 1, 
                       colwave = "black", 
                       coly0 = "lightgrey",
                       colcursor = "red", 
                       bty = "l")
{
    ## ERROR MESSAGES
    if (!is.logical(dB)) stop("'dB' is here a logical.")
    if (slidframe <= 1) stop("'slidframe' is too small")  
    if (slidframe > 90) stop("'slidframe' is too large")
    if (!is.null(envt) && osc) stop("'envt' and 'osc' cannot be used together") 
    ## INPUT
    input <- inputw(wave = wave, f = f, channel = channel) ; wave <- input$w ; f <- input$f ; rm(input)
    ## FROM/TO
    if(!is.null(from)|!is.null(to))
    {
      if(is.null(from) && !is.null(to)) {a<-1; b<-round(to*f)}
      if(!is.null(from) && is.null(to)) {a<-round(from*f)+1; b<-length(wave)}
      if(!is.null(from) && !is.null(to))
        {
          if(from > to) stop("'from' cannot be superior to 'to'")
          if(from==0) {a<-1} else {a<-round(from*f)+1}
          b<-round(to*f)
        }
      wave<-as.matrix(wave[a:b,])
    }
    ## STDFT
    n <- nrow(wave)
    step <- seq(1, n - wl, wl - (ovlp * wl/100))
    lstep <- length(step)
    z <- stdft(wave = wave, f = f, wl = wl, zp = zp, step = step, wn = wn, fftw = fftw)
    if(dB) {z <- 20*log10(z)}  
    y <- seq(0, (f/2) - (f/(wl + zp)), by = f/(wl + zp))/1000
    ## OSCILLO
    if (osc) {
        max <- max(abs(wave))
        losc <- oscillo(wave = wave, f = f, plot=FALSE)
        losc <- cbind(losc, (1:n)/f)
        losc2 <- losc[seq(1, nrow(losc), by=50),]
        rm(losc)
    }
    ## PLOT
    if (plot) {
        pos <- 1:(lstep - floor(lstep * slidframe / 100))
        poslabel <- numeric(lstep)
        poslabel <- (step + ((step[2] - step[1])/2))/f
        z <- t(z)
        ## PLOT FUNCTION
        plot.spec <- function(panel) {
            with(panel, {
                ## LAYOUT
                par(bg=colbg, col=colline, col.axis=colaxis, col.lab=collab,
                    cex.lab=cexlab, font.lab=fontlab, las=1)
                if (osc | !is.null(envt)) {
                    layout(c(1, 2), heights = c(3.5, 1))
                    par(mar = c(4.5, 4, 3, 2))
                }
                ## SPECTRO
                if(is.null(collevels)) {                                     
                    if(dB) {collevels <- seq(-30,0,1)}                        
                    else {collevels <- seq(min(z),max(z),length.out=30)}     
                    }                                                        
       		image(x=poslabel, y=y, z=z, useRaster=TRUE,
                      zlim=c(min(collevels), max(collevels)),                
                      xlab=tlab, ylab=flab,                                  
                      col=palette(length(collevels)),                        
                      xlim=c(poslabel[pos], poslabel[pos + floor(lstep * slidframe / 100) - 1]))
                ## TITLE
                if (title) {
                    nr <- nrow(z)
                    title(main = paste("Position along the signal = ", 
                                       round(poslabel[pos],2), "-",
                                       round(poslabel[pos + floor(lstep * slidframe / 100) - 1],2) ,
                                       "s", sep = " "),
                          col.main = coltitle)
                }
                if (title == FALSE) {title(main = "")}
                if (is.character(title)) {title(main = paste(title), col.main = coltitle)}
                ## OSCILLO / ENVELOPE
                if (osc) {
                    par(mar = c(4.5, 4, 0.5, 2))
                    max <- max(abs(wave))
                    plot(losc2[,2], losc2[,1], type="l", col=colwave, xlab=tlab, ylab=alab, yaxt="n", xaxs="i")
                    abline(h=0, col=coly0)
                    abline(v = c(poslabel[pos], poslabel[pos + floor(lstep * slidframe / 100) - 1]), col = colcursor)
                }
                else if (!is.null(envt)) {
                    par(mar = c(4.5, 4, 0.5, 2))
                    env(wave = wave, f = f, k = 1, j = 1, envt = envt, 
                        msmooth = msmooth, ksmooth = ksmooth, colwave = colwave, 
                        collab = collab, cexlab = cexlab, fontlab = fontlab, 
                        colline = colline, colaxis = colaxis, coly0 = coly0, 
                        bty = bty)
                    abline(v = poslabel[c(pos, (pos+wl))], col = colcursor)
                }
            })
            panel
        }
        ## RPANEL 
        spec.panel <- rpanel::rp.control("Position", size=c(300, 100))
        rpanel::rp.slider(spec.panel, pos, from = 1, to = lstep - floor(lstep * slidframe / 100) - 1, 
                  resolution = 1, title = "Position along the signal", 
                  action = plot.spec)
    }
    ## OUTPUT
    else return(list(time = seq(0, n/f, length.out = length(step)), freq = y, amp = z))
}

################################################################################
##                               ECHO
################################################################################

echo <- function(
                wave,
                f,
                channel = 1,
                amp,
                delay,
                plot= FALSE,
                listen = FALSE,               
                output = "matrix",
                ...
                )

{
  input <- inputw(wave=wave,f=f,channel=channel) ; wave <- input$w ; f <- input$f ; bit <- input$bit ; rm(input)

  wave <- wave/max(abs(wave))
  n <- nrow(wave)
  namp <- length(amp)       # number of amplitude values
  ndelay <- length(delay)   # number of delay values
  delayp <- delay*f         # position of the delays
  delaypn <- delayp[namp]+n # number of samples added by the delays
  
  if(namp!=ndelay) stop("'namp' and 'ndelay' arguments should have the same length")
  if(any(namp)>1)  stop("'namp' argument cannot be > 1")
  if(any(namp)<0)  stop("'namp' argument cannot be negative")

  pulse <- c(1,numeric(delaypn+n-1))
  for(i in 1:namp) pulse[delayp[i]]<-amp[i]
  pulse <- rev(pulse)

  ## test on wave length (n) and pulse length
  ## for a fast convolution we must have the condition: n + pulse length - 1 = x^2 
  ## http://r.789695.n4.nabble.com/stats-convolve-documentation-enhancement-td4670174.html
    fl <- n + length(pulse) - 1
    log.fl <- log2(fl)
    ## if log.fl is not a whole number, that is if the power of 2 condition is false
    ## we have to pad zeros to wave to reach the condition fl = x^2
    if(log.fl != floor(log.fl)){
        ## find the closest above power of 2
        p <- ceiling(log.fl)
        ## number of zeros to padd to wave
        nzp <- 2^p - fl
        wave <- c(wave, rep(0,nzp))
    }

  wave2 <- convolve(wave, pulse, type="open")
  wave2 <- wave2[1:delaypn]
  wave2 <- wave2 - mean(wave2)
  wave2 <- outputw(wave=wave2, f=f, format=output)

  if(plot)
    {
      oscillo(wave=wave2,f=f,...)
      if(listen) {listen(wave2,f=f)}
      invisible(wave2)
    }
  else
    {
      if(listen) {listen(wave2,f=f)}
      return(wave2)
    }
}


################################################################################
##                               ENV                                        
################################################################################

env <- function(wave,
                f,
                channel = 1,
                envt = "hil",
                msmooth = NULL,   
                ksmooth = NULL,
                ssmooth = NULL,
                asmooth = NULL,
                fftw = FALSE,
                norm = FALSE,	
                plot = TRUE,
                k=1,
                j=1,
                ...)

{
    ## STOP MESSAGES
    if(!is.null(msmooth) & !is.null(ksmooth) & !is.null(ssmooth) & !is.null(asmooth)) stop("Please use one of the smoothing arguments, not all of them together.") 
    if(!is.null(msmooth) & !is.null(ksmooth) & !is.null(ssmooth)) stop("Please use one of the smoothing arguments, not three of them.") 
    if(!is.null(msmooth) & !is.null(ssmooth) & !is.null(asmooth)) stop("Please use one of the smoothing arguments, not three of them.") 
    if(!is.null(msmooth) & !is.null(ksmooth) & !is.null(asmooth)) stop("Please use one of the smoothing arguments, not three of them.") 
    if(!is.null(ksmooth) & !is.null(ssmooth) & !is.null(asmooth)) stop("Please use one of the smoothing arguments, not three of them.") 
    if(!is.null(msmooth) & !is.null(ksmooth)) stop("Please use one of the smoothing arguments, not two of them.") 
    if(!is.null(msmooth) & !is.null(ssmooth)) stop("Please use one of the smoothing arguments, not two of them.") 
    if(!is.null(msmooth) & !is.null(asmooth)) stop("Please use one of the smoothing arguments, not two of them.") 
    if(!is.null(ksmooth) & !is.null(ssmooth)) stop("Please use one of the smoothing arguments, not two of them.") 
    if(!is.null(ksmooth) & !is.null(asmooth)) stop("Please use one of the smoothing arguments, not two of them.") 
    if(!is.null(ssmooth) & !is.null(asmooth)) stop("Please use one of the smoothing arguments, not two of them.") 

    ## INPUT
    input <- inputw(wave=wave,f=f,channel=channel) ; wave <- input$w ; f <- input$f ; rm(input)
    n <- nrow(wave)

    if(envt=="hil"){
        wave <- Mod(hilbert(wave,f=f,fftw=fftw))
    } 
    else if(envt=="abs"){
        if(fftw==TRUE && is.null(ksmooth))