--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/nmr_preprocessing/NmrPreprocessing_script.R	Mon Jul 30 10:33:03 2018 -0400
@@ -0,0 +1,1244 @@
+## ==========================
+# Internal functions
+## ==========================
+
+# beginTreatment 
+beginTreatment <- function(name, Signal_data = NULL, Signal_info = NULL, 
+                           force.real = FALSE) {
+  
+  cat("Begin", name, "\n")
+  
+  
+  # Formatting the Signal_data and Signal_info -----------------------
+  
+  vec <- is.vector(Signal_data)
+  if (vec) {
+    Signal_data <- vec2mat(Signal_data)
+  }
+  if (is.vector(Signal_info)) {
+    Signal_info <- vec2mat(Signal_info)
+  }
+  if (!is.null(Signal_data)) {
+    if (!is.matrix(Signal_data)) {
+      stop("Signal_data is not a matrix.")
+    }
+    if (!is.complex(Signal_data) && !is.numeric(Signal_data)) {
+      stop("Signal_data contains non-numerical values.")
+    }
+  }
+  if (!is.null(Signal_info) && !is.matrix(Signal_info)) {
+    stop("Signal_info is not a matrix.")
+  }
+  
+  
+  Original_data <- Signal_data
+  
+  # Extract the real part of the spectrum ---------------------------
+  
+  if (force.real) {
+    if (is.complex(Signal_data)) {
+      Signal_data <- Re(Signal_data)
+    } else {
+      # The signal is numeric Im(Signal_data) is zero anyway so let's avoid
+      # using complex(real=...,imaginary=0) which would give a complex signal
+      # in endTreatment()
+      force.real <- FALSE
+    }
+  }
+  
+  
+  # Return the formatted data and metadata entries --------------------
+  
+  return(list(start = proc.time(), vec = vec, force.real = force.real, 
+              Original_data = Original_data, Signal_data = Signal_data, Signal_info = Signal_info))
+}
+
+# endTreatment 
+endTreatment <- function(name, begin_info, Signal_data) {
+  
+  # begin_info: object outputted from beginTreatment
+  
+  
+  # Formatting the entries and printing process time -----------------------
+  end_time <- proc.time()  # record it as soon as possible
+  start_time <- begin_info[["start"]]
+  delta_time <- end_time - start_time
+  delta <- delta_time[]
+  cat("End", name, "\n")
+  cat("It lasted", round(delta["user.self"], 3), "s user time,", round(delta["sys.self"],3),
+      "s system time and", round(delta["elapsed"], 3), "s elapsed time.\n")
+  
+  
+  if (begin_info[["force.real"]]) {
+    # The imaginary part is left untouched
+    i <- complex(real = 0, imaginary = 1)
+    Signal_data <- Signal_data + i * Im(begin_info[["Original_data"]])
+  }
+  
+  if (begin_info[["vec"]]) {
+    Signal_data <- Signal_data[1, ]
+  }
+  
+  # Return the formatted data and metadata entries --------------------
+  return(Signal_data)
+}
+
+# checkArg 
+checkArg <- function(arg, checks, can.be.null=FALSE) {
+  check.list <- list(bool=c(is.logical, "a boolean"),
+                     int =c(function(x){x%%1==0}, "an integer"),
+                     num =c(is.numeric, "a numeric"),
+                     str =c(is.character, "a string"),
+                     pos =c(function(x){x>0}, "positive"),
+                     pos0=c(function(x){x>=0}, "positive or zero"),
+                     l1 =c(function(x){length(x)==1}, "of length 1")
+  )
+  if (is.null(arg)) {
+    if (!can.be.null) {
+      stop(deparse(substitute(arg)), " is null.")
+    }
+  } else {
+    if (is.matrix(arg)) {
+      stop(deparse(substitute(arg)), " is not scalar.")
+    }
+    for (c in checks) {
+      if (!check.list[[c]][[1]](arg)) {
+        stop(deparse(substitute(arg)), " is not ", check.list[[c]][[2]], ".")
+      }
+    }
+  }
+}
+
+# getArg 
+getArg <- function(arg, info, argname, can.be.absent=FALSE) {
+  if (is.null(arg)) {
+    start <- paste("impossible to get argument", argname, "it was not given directly and");
+    if (!is.matrix(info)) {
+      stop(paste(start, "the info matrix was not given"))
+    }
+    if (!(argname %in% colnames(info))) {
+      if (can.be.absent) {
+        return(NULL)
+      } else {
+        stop(paste(start, "is not in the info matrix"))
+      }
+    }
+    if (nrow(info) < 1) {
+      stop(paste(start, "the info matrix has no row"))
+    }
+    arg <- info[1,argname]
+    if (is.na(arg)) {
+      stop(paste(start, "it is NA in the info matrix"))
+    }
+  }
+  return(arg)
+}
+
+# binarySearch
+binarySearch <- function(a, target, lower = TRUE) {
+  # search the index i in a such that a[i] == target 
+  # if it doesn't exists and lower, it searches the closer a[i] such that a[i] < target
+  # if !lower, it seraches the closer a[i] such that a[i] > target 
+  # a should be monotone but can be increasing or decreasing
+  
+  # if a is increasing INVARIANT: a[amin] < target < a[amax]
+  N <- length(a)
+  if ((a[N] - target) * (a[N] - a[1]) <= 0) {
+    return(N)
+  }
+  if ((a[1] - target) * (a[N] - a[1]) >= 0) {
+    return(1)
+  }
+  amin <- 1
+  amax <- N
+  while (amin + 1 < amax) {
+    amid <- floor((amin + amax)/2)
+    if ((a[amid] - target) * (a[amax] - a[amid]) < 0) {
+      amin <- amid
+    } else if ((a[amid] - target) * (a[amax] - a[amid]) > 0) {
+      amax <- amid
+    } else {
+      # a[amid] == a[amax] or a[amid] == target In both cases, a[amid] ==
+      # target
+      return(amid)
+    }
+  }
+  if (xor(lower, a[amin] > a[amax])) {
+    # (lower && a[amin] < a[amax]) || (!lower && a[min] > a[max]) 
+    # If increasing and we want the lower, we take amin 
+    # If decreasing and we want the bigger, we take amin too
+    return(amin)
+  } else {
+    return(amax)
+  }
+}
+
+# Interpol
+Interpol <- function(t, y) {
+  # y: sample
+  # t : warping function
+  
+  m <- length(y)
+  # t <= m-1
+  # because if t > m-1, y[ti+1] will be NA when we compute g
+  valid <- 1 <= t & t <= m-1 # FIXME it was '<' in Bubble v2
+  s <- (1:m)[valid]
+  ti <- floor(t[s])
+  tr <- t[s] - ti
+  g <- y[ti + 1] - y[ti]
+  f <- y[ti] + tr * g
+  list(f=f, s=s, g=g)
+}
+
+# vec2mat
+vec2mat <- function(vec) {
+  return(matrix(vec, nrow = 1, dimnames = list(c(1), names(vec)))) 
+  
+}
+
+# binarySearch
+binarySearch <- function(a, target, lower = TRUE) {
+  # search the index i in a such that a[i] == target 
+  # if it doesn't exists and lower, it searches the closer a[i] such that a[i] < target
+  # if !lower, it seraches the closer a[i] such that a[i] > target 
+  # a should be monotone but can be increasing or decreasing
+  
+  # if a is increasing INVARIANT: a[amin] < target < a[amax]
+  N <- length(a)
+  if ((a[N] - target) * (a[N] - a[1]) <= 0) {
+    return(N)
+  }
+  if ((a[1] - target) * (a[N] - a[1]) >= 0) {
+    return(1)
+  }
+  amin <- 1
+  amax <- N
+  while (amin + 1 < amax) {
+    amid <- floor((amin + amax)/2)
+    if ((a[amid] - target) * (a[amax] - a[amid]) < 0) {
+      amin <- amid
+    } else if ((a[amid] - target) * (a[amax] - a[amid]) > 0) {
+      amax <- amid
+    } else {
+      # a[amid] == a[amax] or a[amid] == target In both cases, a[amid] ==
+      # target
+      return(amid)
+    }
+  }
+  if (xor(lower, a[amin] > a[amax])) {
+    # (lower && a[amin] < a[amax]) || (!lower && a[min] > a[max]) 
+    # If increasing and we want the lower, we take amin 
+    # If decreasing and we want the bigger, we take amin too
+    return(amin)
+  } else {
+    return(amax)
+  }
+}
+
+
+# indexInterval
+indexInterval <- function (a, from, to, inclusive=TRUE) {
+  # If inclusive and from <= to, we need to take the lower
+  # If not inclusive and from > to, we need to take the lower too
+  lowerFrom <- xor(inclusive, from > to)
+  fromIndex <- binarySearch(a, from, lowerFrom)
+  toIndex <- binarySearch(a, to, !lowerFrom)
+  return(fromIndex:toIndex)
+}
+
+
+
+## ==========================
+# GroupDelayCorrection 
+## ==========================
+GroupDelayCorrection <- function(Fid_data, Fid_info = NULL, group_delay = NULL) {
+  
+    
+    # Data initialisation and checks ----------------------------------------------
+    
+    begin_info <- beginTreatment("GroupDelayCorrection", Fid_data, Fid_info)
+    Fid_data <- begin_info[["Signal_data"]]
+    dimension_names <- dimnames(Fid_data)
+    Fid_info <- begin_info[["Signal_info"]]
+    checkArg(group_delay, c("num", "pos0"), can.be.null = TRUE)
+    # if Fid_info and group_delay are NULL, getArg will generate an error
+    
+    group_delay <- getArg(group_delay, Fid_info, "GRPDLY", can.be.absent = TRUE)
+    
+    if (is.null(group_delay)) {
+      
+      # See DetermineBrukerDigitalFilter.m in matNMR MATLAB library
+      group_delay_matrix <- matrix(c(44.75, 46, 46.311, 33.5, 36.5, 36.53, 66.625, 
+                                     48, 47.87, 59.0833, 50.1667, 50.229, 68.5625, 53.25, 53.289, 60.375, 
+                                     69.5, 69.551, 69.5313, 72.25, 71.6, 61.0208, 70.1667, 70.184, 70.0156, 
+                                     72.75, 72.138, 61.3438, 70.5, 70.528, 70.2578, 73, 72.348, 61.5052, 70.6667, 
+                                     70.7, 70.3789, 72.5, 72.524, 61.5859, 71.3333, NA, 70.4395, 72.25, NA, 
+                                     61.6263, 71.6667, NA, 70.4697, 72.125, NA, 61.6465, 71.8333, NA, 70.4849, 
+                                     72.0625, NA, 61.6566, 71.9167, NA, 70.4924, 72.0313, NA), nrow = 21, 
+                                   ncol = 3, byrow = TRUE, dimnames = list(c(2, 3, 4, 6, 8, 12, 16, 24, 
+                                                                             32, 48, 64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536, 2048), 
+                                                                           c(10, 11, 12)))
+      decim <- Fid_info[1, "DECIM"]
+      dspfvs <- Fid_info[1, "DSPFVS"]
+      if (!(toString(decim) %in% rownames(group_delay_matrix)))  {
+        stop(paste("Invalid DECIM", decim, "it should be one of", rownames(group_delay_matrix)))
+      }
+      if (!(toString(dspfvs) %in% colnames(group_delay_matrix)))  {
+        stop(paste("Invalid DSPFVS", dspfvs, "it should be one of", colnames(group_delay_matrix)))
+      }
+      group_delay <- group_delay_matrix[toString(decim), toString(dspfvs)]
+      if (is.na(group_delay))  {
+        stop(paste("Invalid DECIM", decim, "for DSPFVS", dspfvs))
+      }
+    }
+    m <- ncol(Fid_data)
+    n <- nrow(Fid_data) 
+    
+    # GroupDelayCorrection ----------------------------------------------
+    
+    # We do the shifting in the Fourier domain because the shift can be non-integer.
+    # That way we automatically have the circular behaviour of the shift and the
+    # interpolation if it is non-integer.
+    
+    Spectrum <- t(stats::mvfft(t(Fid_data)))
+    
+    # Spectrum <- FourierTransform(Fid_data, Fid_info)
+    p <- ceiling(m/2)
+    new_index <- c((p + 1):m, 1:p)
+    Spectrum <- Spectrum[,new_index]
+    Spectrum <- matrix(data = Spectrum, ncol = m, nrow = n)
+    
+    Omega <- (0:(m - 1))/m
+    i <- complex(real = 0, imaginary = 1)
+    
+    if (n>1) {
+      Spectrum <- sweep(Spectrum, MARGIN = 2, exp(i * group_delay * 2 * pi * Omega), `*`)
+      Spectrum <- Spectrum[,new_index]
+    }else {
+      Spectrum <- Spectrum* exp(i * group_delay * 2 * pi * Omega)
+      Spectrum <- Spectrum[new_index]
+      Spectrum <- matrix(data = Spectrum, ncol = m, nrow = n)
+    }
+    
+    
+    Fid_data <- t(stats::mvfft(t(Spectrum), inverse = TRUE))/m
+    colnames(Fid_data) <- dimension_names[[2]]
+    rownames(Fid_data) <- dimension_names[[1]]
+    
+    # Data finalisation ----------------------------------------------
+    
+    return(endTreatment("GroupDelayCorrection", begin_info, Fid_data))
+  }
+
+## ==========================
+# SolventSuppression 
+## ==========================
+SolventSuppression <- function(Fid_data, lambda.ss = 1e+06, ptw.ss = TRUE, 
+                               plotSolvent = F, returnSolvent = F) {
+  
+  # Data initialisation and checks ----------------------------------------------
+  
+  begin_info <- beginTreatment("SolventSuppression", Fid_data)
+  Fid_data <- begin_info[["Signal_data"]]
+  checkArg(ptw.ss, c("bool"))
+  checkArg(lambda.ss, c("num", "pos0"))
+  
+  
+  # difsm function definition for the smoother -----------------------------------
+  
+  if (ptw.ss) {
+    # Use of the function in ptw that smoothes signals with a finite difference
+    # penalty of order 2
+    difsm <- ptw::difsm
+  } else  {
+    # Or manual implementation based on sparse matrices for large data series (cf.
+    # Eilers, 2003. 'A perfect smoother')
+    difsm <- function(y, d = 2, lambda)  {
+      
+      m <- length(y)
+      # Sparse identity matrix m x m
+      E <- Matrix::Diagonal(m)
+      D <- Matrix::diff(E, differences = d)
+      A <- E + lambda.ss * Matrix::t(D) %*% D
+      # base::chol does not take into account that A is sparse and is extremely slow
+      C <- Matrix::chol(A)
+      x <- Matrix::solve(C, Matrix::solve(Matrix::t(C), y))
+      return(as.numeric(x))
+    }
+  }
+  
+  # Solvent Suppression ----------------------------------------------
+  
+  n <- dim(Fid_data)[1]
+  if (returnSolvent)  {
+    SolventRe <- Fid_data
+    SolventIm <- Fid_data
+  }
+  for (i in 1:n) {
+    FidRe <- Re(Fid_data[i, ])
+    FidIm <- Im(Fid_data[i, ])
+    solventRe <- difsm(y = FidRe, lambda = lambda.ss)
+    solventIm <- difsm(y = FidIm, lambda = lambda.ss)
+    
+    if (plotSolvent)  {
+      m <- length(FidRe)
+      graphics::plot(1:m, FidRe, type = "l", col = "red")
+      graphics::lines(1:m, solventRe, type = "l", col = "blue")
+      graphics::plot(1:m, FidIm, type = "l", col = "red")
+      graphics::lines(1:m, solventIm, type = "l", col = "blue")
+    }
+    FidRe <- FidRe - solventRe
+    FidIm <- FidIm - solventIm
+    Fid_data[i, ] <- complex(real = FidRe, imaginary = FidIm)
+    if (returnSolvent) {
+      SolventRe[i, ] <- solventRe
+      SolventIm[i, ] <- solventIm
+    }
+  }
+  
+  
+  # Data finalisation ----------------------------------------------
+  
+  Fid_data <- endTreatment("SolventSuppression", begin_info, Fid_data)
+  if (returnSolvent) {
+    return(list(Fid_data = Fid_data, SolventRe = SolventRe, SolventIm = SolventIm))
+  } else {
+    return(Fid_data)
+  }
+}
+
+
+## ==========================
+# Apodization 
+# =============================
+Apodization <- function(Fid_data, Fid_info = NULL, DT = NULL, 
+                        type.apod = c("exp","cos2", "blockexp", "blockcos2", 
+                                      "gauss", "hanning", "hamming"), phase = 0, rectRatio = 1/2, 
+                        gaussLB = 1, expLB = 1, plotWindow = F, returnFactor = F) {
+  
+  # Data initialisation and checks ----------------------------------------------
+  begin_info <- beginTreatment("Apodization", Fid_data, Fid_info)
+  Fid_data <- begin_info[["Signal_data"]]
+  Fid_info <- begin_info[["Signal_info"]]
+  # Data check
+  type.apod <- match.arg(type.apod)
+  checkArg(DT, c("num", "pos"), can.be.null = TRUE)
+  checkArg(phase, c("num"))
+  
+  # Apodization ----------------------------------------------
+  DT <- getArg(DT, Fid_info, "DT")  # Dwell Time
+  m <- ncol(Fid_data)
+  t <- (1:m) * DT  # Time
+  rectSize <- ceiling(rectRatio * m)
+  gaussLB <- (gaussLB/(sqrt(8 * log(2))))
+  # Define the types of apodization:
+  switch(type.apod, exp = {
+    # exponential
+    Factor <- exp(-expLB * t)
+  }, cos2 = {
+    # cos^2
+    c <- cos((1:m) * pi/(2 * m) - phase * pi/2)
+    Factor <- c * c
+  }, blockexp = {
+    # block and exponential
+    Factor <- c(rep.int(1, rectSize), rep.int(0, m - rectSize))
+    # | rectSize | 1 ___________ | \ 0 \____
+    Factor[(rectSize + 1):m] <- exp(-expLB * t[1:(m - rectSize)])
+  }, blockcos2 = {
+    # block and cos^2
+    Factor <- c(rep.int(1, rectSize), rep.int(0, m - rectSize))
+    c <- cos((1:(m - rectSize)) * pi/(2 * (m - rectSize)))
+    Factor[(rectSize + 1):m] <- c * c
+  }, gauss = {
+    # gaussian
+    Factor <- exp(-(gaussLB * t)^2/2)
+    Factor <- Factor/max(Factor)
+  }, hanning = {
+    # Hanning
+    Factor <- 0.5 + 0.5 * cos((1:m) * pi/m - phase * pi)
+  }, hamming = {
+    # Hamming
+    Factor <- 0.54 + 0.46 * cos((1:m) * pi/m - phase * pi)
+  })
+  if (plotWindow) {
+    graphics::plot(1:m, Factor, "l")
+    # dev.off() # device independent, it is the responsability of the
+    # caller to do it
+  }
+  # Apply the apodization factor on the spectra
+  Fid_data <- sweep(Fid_data, MARGIN = 2, Factor, `*`)
+  
+  # Data finalisation ----------------------------------------------
+  Fid_data <- endTreatment("Apodization", begin_info, Fid_data)
+  if (returnFactor) {
+    return(list(Fid_data = Fid_data, Factor = Factor))
+  } else {
+    return(Fid_data)
+  }
+}
+
+
+## ====================================================
+# FourierTransform           
+## ====================================================
+
+
+# fftshift1D2D 
+fftshift1D2D <- function(x) {
+  vec <- F
+  if (is.vector(x)) {
+    x <- vec2mat(x)
+    vec <- T
+  }
+  m <- dim(x)[2]
+  p <- ceiling(m/2)
+  new_index <- c((p + 1):m, 1:p)
+  y <- x[, new_index, drop = vec]
+}
+
+# FourierTransform
+FourierTransform <- function(Fid_data, Fid_info = NULL, SW_h = NULL, SW = NULL, O1 = NULL, reverse.axis = TRUE) {
+  
+  # Data initialisation and checks ----------------------------------------------
+  begin_info <- beginTreatment("FourierTransform", Fid_data, Fid_info)
+  Fid_data <- begin_info[["Signal_data"]]
+  Fid_info <- begin_info[["Signal_info"]]
+  
+  m <- ncol(Fid_data)
+  n <- nrow(Fid_data)
+  
+  if (is.null(SW_h)) {
+    SW_h <- getArg(SW_h, Fid_info, "SW_h")
+  }
+  
+  if (is.null(SW)) {
+    SW <- getArg(SW, Fid_info, "SW")  # Sweep Width in ppm (semi frequency scale in ppm)
+  }
+  
+  
+  if (is.null(O1)) {
+    O1 <- getArg(O1, Fid_info, "O1")
+  }
+  
+  
+  checkArg(reverse.axis, c("bool"))
+  
+  # Fourier Transformation ----------------------------------------------
+  # mvfft does the unnormalized fourier transform (see ?mvfft), so we need divide
+  # by m.  It does not matter a lot in our case since the spectrum will be
+  # normalized.
+  
+  # FT
+  RawSpect_data <- fftshift1D2D(t(stats::mvfft(t(Fid_data))))
+  # recover the frequencies values
+  f <- ((0:(m - 1)) - floor(m/2)) * Fid_info[1, "SW_h"]/(m-1)
+  
+  if(reverse.axis == TRUE) {
+    revind <- rev(1:m)
+    RawSpect_data <- RawSpect_data[,revind] # reverse the spectrum
+  }
+  
+  RawSpect_data <- matrix(RawSpect_data, nrow = n, ncol = m)
+  colnames(RawSpect_data) <- f
+  rownames(RawSpect_data) <- rownames(Fid_data)
+  
+  # PPM conversion ----------------------------------------------
+  
+  # The Sweep Width has to be the same since the column names are the same
+  
+  ppmInterval <- SW/(m-1)
+  
+  O1index = round((m+1)/2+O1*(m - 1) / SW_h)
+  
+  end <- O1index - m
+  start <- O1index -1
+  ppmScale <- (start:end) * ppmInterval
+  RawSpect_data <- matrix(RawSpect_data, nrow = n, ncol =  -(end - start) + 1, dimnames = 
+                            list(rownames(RawSpect_data), ppmScale))
+  
+  
+  # Data finalisation ----------------------------------------------
+  return(endTreatment("FourierTransform", begin_info, RawSpect_data))
+}
+
+## ====================================================
+#   InternalReferencing       
+## ====================================================
+
+InternalReferencing <- function(Spectrum_data, Fid_info, method = c("max", "thres"), 
+                                range = c("nearvalue", "all", "window"), ppm.value = 0, 
+                                direction = "left", shiftHandling = c("zerofilling", "cut", 
+                                                                      "NAfilling", "circular"), c = 2, pc = 0.02, fromto.RC = NULL,
+                                ppm.ir = TRUE, rowindex_graph = NULL) {
+  
+  
+  
+  # Data initialisation and checks ----------------------------------------------
+  
+  begin_info <- beginTreatment("InternalReferencing", Spectrum_data, Fid_info)
+  Spectrum_data <- begin_info[["Signal_data"]]
+  Fid_info <- begin_info[["Signal_info"]]
+  
+  
+  ######## Check input arguments
+  
+  range <- match.arg(range)
+  shiftHandling <- match.arg(shiftHandling)
+  method <- match.arg(method)
+  plots <- NULL
+  
+  checkArg(ppm.ir, c("bool"))
+  checkArg(unlist(fromto.RC), c("num"), can.be.null = TRUE)
+  checkArg(pc, c("num"))
+  checkArg(ppm.value, c("num"))
+  checkArg(rowindex_graph, "num", can.be.null = TRUE)
+  
+  # fromto.RC : if range == "window", 
+  # fromto.RC defines the spectral window where to search for the peak
+  if (!is.null(fromto.RC)) {
+    diff <- diff(unlist(fromto.RC))[1:length(diff(unlist(fromto.RC)))%%2 !=0]
+    for (i in 1:length(diff)) {
+      if (diff[i] >= 0)  {
+        fromto <- c(fromto.RC[[i]][2], fromto.RC[[i]][1])
+        fromto.RC[[i]] <- fromto
+      }
+    }
+  }
+  
+  
+  # findTMSPpeak function ----------------------------------------------
+  # If method == "tresh", findTMSPpeak will find the position of the first 
+  # peak (from left or right) which is higher than a predefined threshold 
+  # and is computed as: c*(cumulated_mean/cumulated_sd)
+  findTMSPpeak <- function(ft, c = 2, direction = "left") {
+    ft <- Re(ft)  # extraction de la partie réelle
+    N <- length(ft)
+    if (direction == "left") {
+      newindex <- rev(1:N)
+      ft <- rev(ft)
+    }
+    thres <- 99999
+    i <- 1000  # Start at point 1000 to find the peak
+    vect <- ft[1:i]
+    
+    while (vect[i] <= (c * thres)) {
+      cumsd <- stats::sd(vect)
+      cummean <- mean(vect)
+      thres <- cummean + 3 * cumsd
+      i <- i + 1
+      vect <- ft[1:i]
+    }
+    if (direction == "left") {
+      v <- newindex[i]
+    } else {v <- i}
+    
+    if (is.na(v))  {
+      warning("No peak found, need to lower the threshold.")
+      return(NA)
+    } else  {
+      # recherche dans les 1% de points suivants du max trouve pour etre au sommet du
+      # pic
+      d <- which.max(ft[v:(v + N * 0.01)])
+      new.peak <- v + d - 1  # nouveau pic du TMSP si d > 0
+      
+      if (names(which.max(ft[v:(v + N * 0.01)])) != names(which.max(ft[v:(v + N * 0.03)])))   {
+        # recherche dans les 3% de points suivants du max trouve pour eviter un faux
+        # positif
+        warning("the TMSP peak might be located further away, increase the threshold to check.")
+      }
+      return(new.peak)
+    }
+  }
+  
+  
+  # Define the search zone  ----------------------------------------
+  
+  n <- nrow(Spectrum_data)
+  m <- ncol(Spectrum_data)
+  
+  # The Sweep Width (SW) has to be the same since the column names are the same
+  SW <- Fid_info[1, "SW"]  # Sweep Width in ppm 
+  ppmInterval <- SW/(m-1)  # size of a ppm interval
+  
+  # range: How the search zone is defined ("all", "nearvalue" or "window")
+  if (range == "all") {
+    
+    Data <- Spectrum_data
+    
+  } else { # range = "nearvalue" or "window"
+    # Need to define colindex (column indexes) to apply indexInterval on it
+    
+    if (range == "nearvalue")  {
+      
+      fromto.RC <- list(c(-(SW * pc)/2 + ppm.value, (SW * pc)/2 + ppm.value))  # automatic fromto values in ppm
+      colindex <- as.numeric(colnames(Spectrum_data))
+      
+    } else {
+      # range == "window"
+      # fromto.RC is already user-defined
+      if (ppm.ir == TRUE)   {
+        colindex <- as.numeric(colnames(Spectrum_data))
+      } else   {
+        colindex <- 1:m
+      }
+    }
+    
+    # index intervals taking into account the different elements in the list fromto.RC
+    Int <- vector("list", length(fromto.RC)) 
+    for (i in 1:length(fromto.RC))  {
+      Int[[i]] <- indexInterval(colindex, from = fromto.RC[[i]][1], 
+                                to = fromto.RC[[i]][2], inclusive = TRUE)
+    }
+    
+    # define Data as the cropped spectrum including the index intervals
+    # outside the research zone, the intensities are set to the minimal 
+    # intensity of the research zone
+    
+    if (n > 1){
+      Data <- apply(Re(Spectrum_data[,unlist(Int)]),1, function(x) rep(min(x), m))
+      Data <- t(Data)
+      Data[,unlist(Int)] <- Re(Spectrum_data[,unlist(Int)])
+    } else {
+      Data <- rep(min(Re(Spectrum_data)) ,m)
+      Data[unlist(Int)] <- Re(Spectrum_data[unlist(Int)])
+    }
+    
+  }
+  
+  
+  # Apply the peak location search method ('thres' or 'max') on spectra
+  # -----------------------------------------------------------------------
+  
+  if (method == "thres") {
+    TMSPpeaks <- apply(Data, 1, findTMSPpeak, c = c, direction = direction)
+  } else { # method == "max
+    TMSPpeaks <- apply(Re(Data), 1, which.max)
+  }
+  
+  
+  # Shift spectra according to the TMSPpeaks found --------------------------------
+  # Depends on the shiftHandling
+  
+  # TMSPpeaks is a column index
+  maxpeak <- max(TMSPpeaks) # max accross spectra
+  minpeak <- min(TMSPpeaks) # min accross spectra
+  
+  
+  if (shiftHandling %in% c("zerofilling", "NAfilling",  "cut")) {
+    fill <- NA
+    if (shiftHandling == "zerofilling")  {
+      fill <- 0
+    }
+    
+    start <-  maxpeak - 1
+    end <- minpeak - m
+    
+    # ppm values of each interval for the whole spectral range of the spectral matrix
+    ppmScale <- (start:end) * ppmInterval
+    
+    # check if ppm.value is in the ppmScale interval
+    if(ppm.value < min(ppmScale) | ppm.value > max(ppmScale)) {
+      warning("ppm.value = ", ppm.value, " is not in the ppm interval [", 
+              round(min(ppmScale),2), ",", round(max(ppmScale),2), "], and is set to its default ppm.value 0")
+      ppm.value = 0
+    }
+    
+    # if ppm.value != 0, ppmScale is adapted
+    ppmScale <- ppmScale + ppm.value
+    
+    # create the spectral matrix with realigned spectra
+    Spectrum_data_calib <- matrix(fill, nrow = n, ncol =  -(end - start) + 1, 
+                                  dimnames = list(rownames(Spectrum_data), ppmScale))
+    
+    # fills in Spectrum_data_calib with shifted spectra
+    for (i in 1:n)  {
+      shift <- (1 - TMSPpeaks[i]) + start
+      Spectrum_data_calib[i, (1 + shift):(m + shift)] <- Spectrum_data[i, ]
+    }
+    
+    if (shiftHandling == "cut")  {
+      Spectrum_data_calib = as.matrix(stats::na.omit(t(Spectrum_data_calib)))
+      Spectrum_data_calib = t(Spectrum_data_calib)
+      base::attr(Spectrum_data_calib, "na.action") <- NULL
+    }
+    
+    
+  } else {
+    # circular
+    start <- 1 - maxpeak
+    end <- m - maxpeak
+    
+    ppmScale <- (start:end) * ppmInterval
+    
+    # check if ppm.value in is the ppmScale interval
+    if(ppm.value < min(ppmScale) | ppm.value > max(ppmScale)) {
+      warning("ppm.value = ", ppm.value, " is not in the ppm interval [", 
+              round(min(ppmScale),2), ",", round(max(ppmScale),2), "], and is set to its default ppm.value 0")
+      ppm.value = 0
+    }
+    
+    # if ppm.value != 0, ppmScale is adapted
+    ppmScale <- ppmScale + ppm.value
+    
+    # create the spectral matrix with realigned spectra
+    Spectrum_data_calib <- matrix(nrow=n, ncol=end-start+1,
+                                  dimnames=list(rownames(Spectrum_data), ppmScale))
+    
+    # fills in Spectrum_data_calib with shifted spectra
+    for (i in 1:n) {
+      shift <- (maxpeak-TMSPpeaks[i])
+      Spectrum_data_calib[i,(1+shift):m] <- Spectrum_data[i,1:(m-shift)]
+      if (shift > 0) {
+        Spectrum_data_calib[i,1:shift] <- Spectrum_data[i,(m-shift+1):m]
+      }
+    }
+  }
+  
+  
+  
+  
+  # Plot of the spectra (depending on rowindex_graph) ---------------------------------------------------
+  
+  ppm = xstart = value = xend = Legend = NULL # only for R CMD check
+  
+  
+  # with the search zone for TMSP and the location of the peaks just found
+  if (!is.null(rowindex_graph)) {
+    
+    if (range == "window")  {
+      if (ppm.ir == TRUE)   {
+        fromto <- fromto.RC
+      } else  {
+        fromto <- list()
+        idcol <- as.numeric(colnames(Spectrum_data))
+        for (i in 1:length(fromto.RC)) {
+          fromto[[i]] <- as.numeric(colnames(Spectrum_data))[fromto.RC[[i]]]
+        }
+      }
+    } else {
+      fromto <- fromto.RC
+    }
+    
+    # TMSPloc in ppm
+    TMSPloc <- as.numeric(colnames(Spectrum_data))[TMSPpeaks[rowindex_graph]]
+    
+    # num plot per window
+    num.stacked <- 6
+    
+    # rectanglar bands of color for the search zone
+    rects <- data.frame(xstart = sapply(fromto, function(x) x[[1]]), 
+                        xend = sapply(fromto, function(x) x[[2]]), 
+                        Legend = "Peak search zone and location")
+    
+    # vlines for TMSP peak
+    addlines <- data.frame(rowname = rownames(Spectrum_data)[rowindex_graph],TMSPloc)
+    
+    nn <- length(rowindex_graph)
+    i <- 1
+    j <- 1
+    plots <- vector(mode = "list", length = ceiling(nn/num.stacked))
+    
+    while (i <= nn) {
+      
+      last <- min(i + num.stacked - 1, nn)
+      
+      melted <- reshape2::melt(Re(Spectrum_data[i:last, ]), 
+                               varnames = c("rowname", "ppm"))
+      
+      plots[[j]] <- ggplot2::ggplot() + ggplot2::theme_bw() + 
+        ggplot2::geom_line(data = melted, 
+                           ggplot2::aes(x = ppm, y = value)) + 
+        ggplot2::geom_rect(data = rects, ggplot2::aes(xmin = xstart, xmax = xend, 
+                                                      ymin = -Inf, ymax = Inf, fill = Legend), alpha = 0.4) + 
+        ggplot2::facet_grid(rowname ~ ., scales = "free_y") + 
+        ggplot2::theme(legend.position = "none") + 
+        ggplot2::geom_vline(data = addlines, ggplot2::aes(xintercept = TMSPloc), 
+                            color = "red", show.legend = TRUE) + 
+        ggplot2::ggtitle("Peak search zone and location") + 
+        ggplot2::theme(legend.position = "top", legend.text = ggplot2::element_text())
+      
+      
+      
+      if ((melted[1, "ppm"] - melted[(dim(melted)[1]), "ppm"]) > 0) {
+        plots[[j]] <- plots[[j]] + ggplot2::scale_x_reverse()
+      }
+      
+      i <- last + 1
+      j <- j + 1
+    }
+    
+    plots
+  }
+  
+  
+  # Return the results ----------------------------------------------
+  Spectrum_data <- endTreatment("InternalReferencing", begin_info, Spectrum_data_calib)
+  
+  if (is.null(plots)) {
+    return(Spectrum_data)
+  } else {
+    return(list(Spectrum_data = Spectrum_data, plots = plots))
+  }
+  
+}
+
+## ====================================================
+# ZeroOrderPhaseCorrection       
+## ====================================================
+
+ZeroOrderPhaseCorrection <- function(Spectrum_data, type.zopc = c("rms", "manual", "max"), 
+                                     plot_rms = NULL, returnAngle = FALSE, createWindow = TRUE, 
+                                     angle = NULL, plot_spectra = FALSE,  
+                                     ppm.zopc = TRUE, exclude.zopc = list(c(5.1,4.5))) {
+  
+  
+  # Data initialisation and checks ----------------------------------------------
+  
+  # Entry arguments definition:
+  # plot_rms : graph of rms criterion returnAngle : if TRUE, returns avector of
+  # optimal angles createWindow : for plot_rms plots angle : If angle is not NULL,
+  # spectra are rotated according to the angle vector values
+  # plot_spectra : if TRUE, plot rotated spectra  
+  
+  
+  
+  begin_info <- beginTreatment("ZeroOrderPhaseCorrection", Spectrum_data)
+  Spectrum_data <- begin_info[["Signal_data"]]
+  n <- nrow(Spectrum_data)
+  m <- ncol(Spectrum_data)
+  
+  rnames <- rownames(Spectrum_data)
+  
+  # Check input arguments
+  type.zopc <- match.arg(type.zopc)
+  checkArg(ppm.zopc, c("bool"))
+  checkArg(unlist(exclude.zopc), c("num"), can.be.null = TRUE)
+  
+  
+  # type.zopc in c("max", "rms") -----------------------------------------
+  if (type.zopc %in% c("max", "rms")) {
+    # angle is found by optimization
+    
+    # rms function to be optimised
+    rms <- function(ang, y, meth = c("max", "rms"))  {
+      # if (debug_plot) { graphics::abline(v=ang, col='gray60') }
+      roty <- y * exp(complex(real = 0, imaginary = ang))  # spectrum rotation
+      Rey <- Re(roty)
+      
+      if (meth == "rms")  {
+        ReyPos <- Rey[Rey >= 0]  # select positive intensities
+        POSss <- sum((ReyPos)^2, na.rm = TRUE)  # SS for positive intensities
+        ss <- sum((Rey)^2, na.rm = TRUE)  #  SS for all intensities
+        return(POSss/ss)  # criterion : SS for positive values / SS for all intensities 
+      } else  {
+        maxi <- max(Rey, na.rm = TRUE)
+        return(maxi)
+      }
+    }
+    
+    
+    # Define the interval where to search for (by defining Data)
+    if (is.null(exclude.zopc)) {
+      Data <- Spectrum_data
+    } else  {
+      
+      # if ppm.zopc == TRUE, then exclude.zopc is in the colnames values, else, in the column
+      # index
+      if (ppm.zopc == TRUE)  {
+        colindex <- as.numeric(colnames(Spectrum_data))
+      } else  {
+        colindex <- 1:m
+      }
+      
+      # Second check for the argument exclude.zopc
+      diff <- diff(unlist(exclude.zopc))[1:length(diff(unlist(exclude.zopc)))%%2 !=0]
+      for (i in 1:length(diff)) {
+        if (ppm.zopc == TRUE & diff[i] >= 0)  {
+          stop(paste("Invalid region removal because from <= to in ppm.zopc"))
+        } else if (ppm.zopc == FALSE & diff[i] <= 0) {stop(paste("Invalid region removal because from >= to in column index"))}
+      }
+      
+      
+      Int <- vector("list", length(exclude.zopc))
+      for (i in 1:length(exclude.zopc))  {
+        Int[[i]] <- indexInterval(colindex, from = exclude.zopc[[i]][1], 
+                                  to = exclude.zopc[[i]][2], inclusive = TRUE)
+      }
+      
+      vector <- rep(1, m)
+      vector[unlist(Int)] <- 0
+      if (n > 1)  {
+        Data <- sweep(Spectrum_data, MARGIN = 2, FUN = "*", vector)  # Cropped_Spectrum
+      } else   {
+        Data <- Spectrum_data * vector
+      }  # Cropped_Spectrum
+    }
+    
+    
+    # angles computation
+    Angle <- c()
+    for (k in 1:n)
+    {
+      # The function is rms is periodic (period 2pi) and it seems that there is a phase
+      # x such that rms is unimodal (i.e. decreasing then increasing) on the interval
+      # [x; x+2pi].  However, if we do the optimization for example on [x-pi; x+pi],
+      # instead of being decreasing then increasing, it might be increasing then
+      # decreasing in which case optimize, thinking it is a valley will have to choose
+      # between the left or the right of this hill and if it chooses wrong, it will end
+      # up at like x-pi while the minimum is close to x+pi.
+      
+      # Supposing that rms is unimodal, the classical 1D unimodal optimization will
+      # work in either [-pi;pi] or [0;2pi] (this is not easy to be convinced by that I
+      # agree) and we can check which one it is simply by the following trick
+      
+      f0 <- rms(0, Data[k, ],meth = type.zopc)
+      fpi <- rms(pi, Data[k, ], meth = type.zopc)
+      if (f0 < fpi) {
+        interval <- c(-pi, pi)
+      } else {
+        interval <- c(0, 2 * pi)
+      }
+      
+      # graphs of rms criteria
+      debug_plot <- F  # rms should not plot anything now, only when called by optimize
+      if (!is.null(plot_rms) && rnames[k] %in% plot_rms) {
+        x <- seq(min(interval), max(interval), length.out = 100)
+        y <- rep(1, 100)
+        for (K in (1:100))   {
+          y[K] <- rms(x[K], Data[k, ],  meth = type.zopc)
+        }
+        if (createWindow == TRUE)  {
+          grDevices::dev.new(noRStudioGD = FALSE)
+        }
+        graphics::plot(x, y, main = paste("Criterion maximization \n", 
+                                          rownames(Data)[k]), ylim = c(0, 1.1),
+                       ylab = "positiveness criterion", xlab = "angle ")
+        debug_plot <- T
+      }
+      
+      # Best angle
+      best <- stats::optimize(rms, interval = interval, maximum = TRUE, 
+                              y = Data[k,],  meth = type.zopc)
+      ang <- best[["maximum"]]
+      
+      
+      if (debug_plot)  {
+        graphics::abline(v = ang, col = "black")
+        graphics::text(x = (ang+0.1*ang), y = (y[ang]-0.1*y[ang]), labels = round(ang, 3))
+      }
+      
+      # Spectrum rotation
+      Spectrum_data[k, ] <- Spectrum_data[k, ] * exp(complex(real = 0, imaginary = ang))
+      Angle <- c(Angle, ang)
+    }
+    
+    
+    
+    
+  } else {
+    # type.zopc is "manual" -------------------------------------------------------
+    # if Angle is already specified and no optimisation is needed
+    Angle <- angle
+    
+    if (!is.vector(angle)) {
+      stop("angle is not a vector")
+    }
+    
+    if (!is.numeric(angle))  {
+      stop("angle is not a numeric")
+    }
+    
+    if (length(angle) != n) {
+      stop(paste("angle has length", length(angle), "and there are", n, "spectra to rotate."))
+    }
+    for (k in 1:n)  {
+      Spectrum_data[k, ] <- Spectrum_data[k, ] * exp(complex(real = 0, imaginary = - angle[k]))
+    }
+  }
+  
+  
+  #  Draw spectra
+  if (plot_spectra == TRUE) {
+    nn <- ceiling(n/4)
+    i <- 1
+    for (k in 1:nn)  {
+      if (createWindow == TRUE)  {
+        grDevices::dev.new(noRStudioGD = FALSE)
+      }
+      graphics::par(mfrow = c(4, 2))
+      while (i <= n)   {
+        last <- min(i + 4 - 1, n)
+        graphics::plot(Re(Spectrum_data[i, ]), type = "l", ylab = "intensity", 
+                       xlab = "Index", main = paste0(rownames(Spectrum_data)[i], " - Real part"))
+        graphics::plot(Im(Spectrum_data[i, ]), type = "l", ylab = "intensity", 
+                       xlab = "Index", main = paste0(rownames(Spectrum_data)[i], " - Imaginary part"))
+        i <- i + 1
+      }
+      i <- last + 1
+    }
+  }
+  
+  
+  # Data finalisation ----------------------------------------------
+  
+  Spectrum_data <- endTreatment("ZeroOrderPhaseCorrection", begin_info, Spectrum_data)
+  if (returnAngle) {
+    return(list(Spectrum_data = Spectrum_data, Angle = Angle))
+  } else {
+    return(Spectrum_data)
+  }
+}
+
+
+## ====================================================
+# Baseline Correction   
+## ====================================================
+BaselineCorrection <- function(Spectrum_data, ptw.bc = TRUE, maxIter = 42, 
+                               lambda.bc = 1e+07, p.bc = 0.05, eps = 1e-08, 
+                               ppm.bc = TRUE, exclude.bc = list(c(5.1,4.5)),
+                               returnBaseline = F) {
+  
+  # Data initialisation ----------------------------------------------
+  begin_info <- beginTreatment("BaselineCorrection", Spectrum_data, force.real = T)
+  Spectrum_data <- begin_info[["Signal_data"]]
+  p <- p.bc
+  lambda <- lambda.bc
+  n <- dim(Spectrum_data)[1]
+  m <- dim(Spectrum_data)[2]
+  
+  
+  # Data check
+  checkArg(ptw.bc, c("bool"))
+  checkArg(maxIter, c("int", "pos"))
+  checkArg(lambda, c("num", "pos0"))
+  checkArg(p.bc, c("num", "pos0"))
+  checkArg(eps, c("num", "pos0"))
+  checkArg(returnBaseline, c("bool"))
+  checkArg(ppm.bc, c("bool"))
+  checkArg(unlist(exclude.bc), c("num"), can.be.null = TRUE)
+  
+  # Define the interval where to search for (by defining Data)
+  if (is.null(exclude.bc)) {
+    exclude_index <- NULL
+  } else  {
+    # if ppm.bc == TRUE, then exclude.bc is in the colnames values, else, in the column
+    # index
+    if (ppm.bc == TRUE)  {
+      colindex <- as.numeric(colnames(Spectrum_data))
+    } else  {
+      colindex <- 1:m
+    }
+    
+    Int <- vector("list", length(exclude.bc))
+    for (i in 1:length(exclude.bc))  {
+      Int[[i]] <- indexInterval(colindex, from = exclude.bc[[i]][1], 
+                                to = exclude.bc[[i]][2], inclusive = TRUE)
+    }
+    exclude_index <- unlist(Int)
+  }
+  
+  # Baseline Correction implementation definition ----------------------
+  
+  # 2 Ways: either use the function asysm from the ptw package or by 
+  # built-in functions 
+  if (ptw.bc) {
+    asysm <- ptw::asysm
+  } else {
+    difsmw <- function(y, lambda, w, d) {
+      # Weighted smoothing with a finite difference penalty cf Eilers, 2003.
+      # (A perfect smoother) 
+      # y: signal to be smoothed 
+      # lambda: smoothing parameter 
+      # w: weights (use0 zeros for missing values) 
+      # d: order of differences in penalty (generally 2)
+      m <- length(y)
+      W <- Matrix::Diagonal(x=w)
+      E <- Matrix::Diagonal(m)
+      D <- Matrix::diff(E, differences = d)
+      C <- Matrix::chol(W + lambda * t(D) %*% D)
+      x <- Matrix::solve(C, Matrix::solve(t(C), w * y))
+      return(as.numeric(x))
+      
+    }
+    asysm <- function(y, lambda, p, eps, exclude_index) {
+      # Baseline estimation with asymmetric least squares
+      # y: signal
+      # lambda: smoothing parameter (generally 1e5 to 1e8)
+      # p: asymmetry parameter (generally 0.001)
+      # d: order of differences in penalty (generally 2)
+      # eps: 1e-8 in ptw package
+      m <- length(y)
+      w <- rep(1, m)
+      i <- 1
+      repeat {
+        z <- difsmw(y, lambda, w, d = 2)
+        w0 <- w
+        p_vect <- rep((1-p), m) # if y <= z + eps
+        p_vect[y > z + eps | y < 0] <- p  # if y > z + eps | y < 0
+        if(!is.null(exclude_index)){
+          p_vect[exclude_index] <- 0 # if exclude area
+        }
+        
+        w <- p_vect  
+        # w <- p * (y > z + eps | y < 0) + (1 - p) * (y <= z + eps)
+        
+        if (sum(abs(w - w0)) == 0) {
+          break
+        }
+        i <- i + 1
+        if (i > maxIter) {
+          warning("cannot find Baseline estimation in asysm")
+          break
+        }
+      }
+      return(z)
+    }
+  }
+  
+  # Baseline estimation ----------------------------------------------
+  Baseline <- matrix(NA, nrow = nrow(Spectrum_data), ncol = ncol(Spectrum_data))
+  
+  # for (k in 1:n) {
+  # Baseline[k, ] <- asysm(y = Spectrum_data[k, ], lambda = lambda, p = p, eps = eps)
+  
+  if (ptw.bc ){
+    Baseline <- apply(Spectrum_data,1, asysm, lambda = lambda, p = p, 
+                      eps = eps)
+  }else {
+    Baseline <- apply(Spectrum_data,1, asysm, lambda = lambda, p = p, 
+                      eps = eps, exclude_index = exclude_index)
+  }
+  
+  
+  Spectrum_data <- Spectrum_data - t(Baseline)
+  # }
+  
+  # Data finalisation ----------------------------------------------
+  Spectrum_data <- endTreatment("BaselineCorrection", begin_info, Spectrum_data)  # FIXME create removeImaginary filter ??
+  
+  if (returnBaseline) {
+    return(list(Spectrum_data = Spectrum_data, Baseline = Baseline))
+  } else {
+    return(Spectrum_data)
+  }
+}
+
+
+
+## ====================================================
+# NegativeValuesZeroing   
+## ====================================================
+
+NegativeValuesZeroing <- function(Spectrum_data) {
+  # Data initialisation and checks ----------------------------------------------
+  begin_info <- beginTreatment("NegativeValuesZeroing", Spectrum_data, force.real = T)
+  Spectrum_data <- begin_info[["Signal_data"]]
+  
+  # NegativeValuesZeroing ----------------------------------------------
+  Spectrum_data[Spectrum_data < 0] <- 0
+  
+  # Data finalisation ----------------------------------------------
+  return(endTreatment("NegativeValuesZeroing", begin_info, Spectrum_data))
+}
+
+