# Computing the Fast Fourier Transform (slightly modified function of the signal R-package)

# INPUT: 	wave: object of Wave-class (see tuneR package)
#			ctrl: control object

# OUTPUT: 	ampl: matrix of signal amplitudes
#			spec: matrix of spectral magnitudes
#			spec.compl: matrix of complex Fourier coefficients
#			phase: matrix of spectral phase
#			f: vector of frequencies for Fourier coefficients
#			win.start.sample: windows starting sample numbers (vector)
#			win.start.time: windows starting times (vector)

fftArray=function (wave, ctrl) 
{
  if (!class(wave)=="Wave") 
    stop("'wave' has to be a 'wave' object from tuneR package.")
  x=wave@left
  win.size = ctrl$window.size
  hop.size = ctrl$hop.size
  samp.rate= ctrl$sampling.rate
  
  #### Normaler weg
  #if (length(x) > win.size) 
  #  offset = seq(1, length(x) - win.size, by = hop.size)
  #else offset = 1
  #S = matrix(0, win.size, length(offset))
  
  #### Boeckers Weg
  frames = length(x)/ hop.size
  seek = (1:frames)*hop.size - win.size
  # Addiere Nullen links vom Signal!
  if (seek[1]<0) {x=c(rep(0,abs(seek[1])),x)}
  offset = seq(1, length(x) - win.size, by = hop.size)
  S = matrix(0, win.size, length(offset))

  
  if (ctrl$window.fun=="rect") window.fun=boxcar
  if (ctrl$window.fun=="hamm") window.fun=hamming
  if (ctrl$window.fun=="hann") window.fun=hanning
  if (ctrl$window.fun=="black") window.fun=blackman
  if (ctrl$window.fun=="gauss") window.fun=gausswin
  

    window.vector=window.fun(win.size)
    
  ## Einschub
  #i= 1 
  #plot(window.vector)
  #plot(x[offset[i]:(offset[i] + win.size - 1)])
  #print(offset[i])
  #print((offset[i] + win.size - 1))
  #plot(x[offset[i]:(offset[i] + win.size - 1)] * window.vector)
  
  for (i in seq_along(offset)) {
    S[1:win.size, i] = x[offset[i]:(offset[i] + win.size - 1)] * window.vector
  }
  A=S
  S <- mvfft(S)
 
  half.win.size=win.size/2
  S.half <- S[1:half.win.size, ]
  f = (1:half.win.size) * samp.rate /  win.size
  win.start.sample = offset
  win.start.time = win.start.sample / samp.rate
  
  
  # spec ist die magnitude (also abs. Koeffiziente) und phase die Phase
  res <- list(ampl=A, spec = abs(S.half), spec.compl = S.half, phase = Arg(S.half), f = f, win.start.sample = win.start.sample, 
              win.start.time =  win.start.time)
  class(res) <- "fftArray"
  return(res)
}


 stand_01=function(v){
   return( (v-min(v)) / (max(v)-min(v)) )
 }


# print methode logarithmiert und auf 0-1 standardisiert!
print.fftArray=function (x, xlab = "Zeit in s", ylab = "Frequenz", leg=TRUE, log=TRUE,...) 
{
  library(fields)
  abs_v=abs(as.vector(x$spec))
  
  # replace 0 with 0.00001
  w=which(abs_v==0)
  abs_v[w]=0.00001
  
  if(log==TRUE)
    SPEC=stand_01(log(abs_v))
  else
    SPEC=stand_01(abs_v)
  if (any(SPEC<0)) 
    stop("Logarithming cases negative values. Try with log=FALSE")
  
  SPEC_M=matrix(SPEC,ncol=ncol(x$spec))
  if (leg==TRUE)
    image.plot(x$win.start.time, x$f, t(SPEC_M), col = topo.colors(20),xlab=xlab,ylab=ylab,...)
  else
    image(x$win.start.time, x$f, t(SPEC_M), col = topo.colors(20),xlab=xlab,ylab=ylab,...)
}
                                       