# computeSELFeatures computes the selected odf features 

# INPUT: 	fa: output of fftArray-function
#			ctrl: ctrl object (a list)
# 			sel.feat: a vector of features with has to be calculated for odfMatrix.

# OUTPUT:	odf.Matrix.sel: odfMatrix of selected features
#			

# sel.feat: nur Namen!
computeSELFeatures=function(fa, ctrl, sel.feat)
{

  S = fa$spec
  k = dim(S)[1] # nr. of frequency bins
  m = dim(S)[2] # nr. of windows
  n = m
  N = k

  NAME = NULL
  # 1
  if (any (sel.feat == "ZCRate")){
    ZCR=ZerroCrossingRate(fa, ctrl)
    ZCRD=abs(diff(c(0,ZCR)))
  }
  
  if (any(sel.feat == "AMax") |  any(sel.feat == "AMaxA")){
    Max.of.wind=apply(fa$ampl, 2, function(x) max(abs(x)))
    AMD=diff(c(0, Max.of.wind))
    AMDA=abs(AMD)
  }
  
  if (any(sel.feat == "AEnerg") |any(sel.feat == "AEnergA") ){
    AE=AmplitEnergy(fa,ctrl)
    AED=diff(c(0, AE)) 
    AEDA=abs(AED) 
  }
  
  if (any(sel.feat == "HFCont") | any(sel.feat == "HFContA")) {
    HFC=HighFreqCont(fa, ctrl)
    HFCD=diff(c(0, HFC)) 
    HFCDA=abs(HFCD)
  }
  
  if (any(sel.feat == "GFCont") | any(sel.feat == "GFContA")) {
    HFCG=HighFreqContGaus(fa, ctrl)
    HFCGD=diff(c(0, HFCG)) 
    HFCGDA=abs(HFCGD)
  }
  


  if (any(sel.feat == "SCentr") | any(sel.feat == "SSpread") | any(sel.feat == "SSkew")) {
    SC=SpecCentroid(fa, ctrl)
    SCD=abs(diff(c(0,SC)))
  
  
 
      SS=numeric(m)
    SC=SpecCentroid(fa, ctrl)

    for (i in 1:m) {
      abs.S=abs(S[,i])
      SS[i]=sqrt( sum(abs.S*((1:k)-SC[i])^2)) /
        sqrt(sum(abs.S))
    }
    SSD=abs(diff(c(0,SS)))
  
  
    SSk=numeric(m)
    SC=SpecCentroid(fa, ctrl)
    for (i in 1:m) {
      abs.S=abs(S[,i])
      SSk[i]=sum(abs.S*((1:k)-SC[i])^3) /
        (SS[i]^3*(sum(abs.S)))
    }
    SSkD=abs(diff(c(0,SSk)))
  }
  
  
  if (any(sel.feat == "SFlux")) {
    SF = SpecFlux(fa, ctrls) 
  } 
  
  if (any(sel.feat == "SFluxEucl")) {
    SFE = SpecFluxEuclid(fa, ctrls) 
  }
  
  
  
  # Phase deviation and Co
  if (any(sel.feat == "PDev") | any(sel.feat == "WPDev") | any(sel.feat == "CDom") | any(sel.feat == "RCDom")) {
    
	S.c = fa$spec.compl
    S.c=cbind(rep(0,dim(S.c)[1]),rep(0,dim(S.c)[1]),S.c)
    
    # intialization
    phi=Arg(S.c)
    phi_2=phi
    PD=numeric(n+2)
    WPD=numeric(n+2)
    S.estim=S.c
    CD=numeric(n+2)
    RCD=numeric(n+2)
    phi.est=phi
    
    for (i in 3:(n+2)) {
      phi_2[,i]=phi[,i]-2*phi[,(i-1)]+phi[,(i-2)]
      phi.est[,i]=2*phi[,(i-1)]-phi[,(i-2)]
      PD[i]=(1/dim(S.c)[1])*sum(abs(phi_2[,i]))
      WPD[i]=(1/dim(S.c)[1])*sum(abs(S.c[,i]*phi_2[,i]))
      S.estim[,i]=complex(modulus = Mod(S.c[,(i-1)]), argument = phi.est[,i])
      CD[i]=(1/dim(S.c)[1])*sum(abs(S.c[,i]-S.estim[,i]))
      
    }
    
    a.S = abs(S.c)-abs(S.estim)
    w.pos = a.S > 0
    A = abs(S.c-S.estim) * w.pos
    RCD = apply(A, 2, sum) * (1/dim(S.c)[1])
    
    
    PD = PD[-(1:2)]
    WPD = WPD[-(1:2)]
    CD = CD[-(1:2)]
    RCD = RCD[-(1:2)]
  }
 
  #names.all= c("ZCRate", "AMax", "AMaxA","AEnerg", "AEnergA", "HFCont", "HFContA", "GFCont", "GFContA", "SCentr", "SSpread", "SSkew", "SFlux",
  #"SFluxEucl", "PDev", "WPDev", "CDom", "RCDom")
  #names.sel = intersect(names.all, sel.feat) # namen der geforderten Vekroren
  
match.ofd = function(v)
{

  s = sapply(v, function (x ) switch(x,  "ZCRate" = ZCRD, "AMax"=AMD, "AMaxA"=AMDA,"AEnerg"=AED, "AEnergA"=AEDA, "HFCont"=HFCD, 
                                     "HFContA"=HFCDA, "GFCont" = HFCGD, "GFContA" = HFCGDA, "SCentr" = SCD, "SSpread"=SSD, "SSkew"=SSkD, "SFlux"=SF,
                                     "SFluxEucl"=SFE, "PDev" = PD, "WPDev" = WPD, "CDom" = CD, "RCDom"= RCD))
  s = t(as.data.frame(s))
 
  return(s)
}

  
 
  odf.Matrix.sel = match.ofd(sel.feat)
  return(odf.Matrix.sel)
}

# not important here
computeALLFeatures=function(fa, ctrl)
{
  NAME = NULL
  # 1
  ZCR=ZerroCrossingRate(fa, ctrl)
  ZCRD=abs(diff(c(0,ZCR)))
  NAME =c(NAME, "ZCRate")

  
  Max.of.wind=apply(fa$ampl, 2, function(x) max(abs(x)))
  AMD=diff(c(0, Max.of.wind))
  NAME =c(NAME, "AMax")
  
  AMDA=abs(AMD)
  NAME =c(NAME, "AMaxA") 
  
  
  AE=AmplitEnergy(fa,ctrl)
  AED=diff(c(0, AE)) 
  NAME =c(NAME, "AEnerg") 
  
  
  AEDA=abs(AED) 
  NAME =c(NAME, "AEnergA") 
  
  
  HFC=HighFreqCont(fa, ctrl)
  HFCD=diff(c(0, HFC)) 
  NAME =c(NAME, "HFCont") 
  
  HFCDA=abs(HFCD)
  NAME =c(NAME, "HFContA") 

  
  HFCG=HighFreqContGaus(fa, ctrl)
  HFCGD=diff(c(0, HFCG)) 
  NAME =c(NAME, "GFCont") 
  
  
  HFCGDA=abs(HFCGD)
  NAME =c(NAME, "GFContA") 
  
  
  # Spectral centroid differenz
  SC=SpecCentroid(fa, ctrl)
  SCD=abs(diff(c(0,SC)))
  NAME =c(NAME, "SCentr") 
  
  
  S = fa$spec
  k = dim(S)[1] # nr. of frequency bins
  m = dim(S)[2] # nr. of windows
  
  SS=numeric(m)
  for (i in 1:m) {
    abs.S=abs(S[,i])
    SS[i]=sqrt( sum(abs.S*((1:k)-SC[i])^2)) /
      sqrt(sum(abs.S))
  }
  SSD=abs(diff(c(0,SS)))
  NAME =c(NAME, "SSpread") 

  
  
  SSk=numeric(m)
  for (i in 1:m) {
    abs.S=abs(S[,i])
    SSk[i]=sum(abs.S*((1:k)-SC[i])^3) /
      (SS[i]^3*(sum(abs.S)))
  }
  SSkD=abs(diff(c(0,SSk)))
  NAME =c(NAME, "SSkew") 
  
  
  
  SF = SpecFlux(fa, ctrls) 
  NAME =c(NAME, "SFlux") 
  
  SFE = SpecFluxEuclid(fa, ctrls) 
  NAME =c(NAME, "SFluxEucl") 
  

  # Phase deviation
  S.c = fa$spec.compl

	
  n = dim(S.c)[1]  # number of frequency bins

  N = dim(S.c)[2] # number of windows

 
  S.c=cbind(rep(0,n),rep(0,n),S.c)
  #print(dim(S.c))
  
  # intialization
  phi=Arg(S.c)
  #print("dim(phi)")
  #print(dim(phi))
  phi_2=phi
  PD=numeric(N+2)
  WPD=numeric(N+2)
  S.estim=S.c
  CD=numeric(N+2)
  RCD=numeric(N+2)
  phi.est=phi
  
  for (i in 3:(N+2)) {
    phi_2[,i]=phi[,i]-2*phi[,(i-1)]+phi[,(i-2)]
    phi.est[,i]=2*phi[,(i-1)]-phi[,(i-2)]
    PD[i]=(1/dim(S.c)[1])*sum(abs(phi_2[,i]))
    WPD[i]=(1/dim(S.c)[1])*sum(abs(S.c[,i]*phi_2[,i]))
    S.estim[,i]=complex(modulus = Mod(S.c[,(i-1)]), argument = phi.est[,i])
    CD[i]=(1/dim(S.c)[1])*sum(abs(S.c[,i]-S.estim[,i]))
    
    
  }

  
  a.S = abs(S.c)-abs(S.estim)
  w.pos = a.S > 0
  A = abs(S.c-S.estim) * w.pos
  RCD = apply(A, 2, sum) * (1/dim(S.c)[1])
  
  
  PD = PD[-(1:2)]
  WPD = WPD[-(1:2)]
  CD = CD[-(1:2)]
  RCD = RCD[-(1:2)]

  

  NAME =c(NAME, "PDev", "WPDev", "CDom", "RCDom") 
  
  odf.Matrix = rbind(ZCRD, AMD, AMDA, AED, AEDA, HFCD, HFCDA, HFCGD, HFCGDA, SCD, SSD, SSkD, SF, SFE, PD, WPD, CD, RCD)
  rownames(odf.Matrix)=NAME
  return(odf.Matrix)
}

odf_function=function (name,fun) 
{
  structure(fun, name = name)
}


ZerroCrossingRate=function(fa,ctrl)
{
  A=fa$ampl
  A_mov=rbind(A[-1,],rep(0,ncol(A)))
  AA=A*A_mov
  ZCR=apply(AA,2,function(x)sum(x<0))
  return(ZCR)
}


ZerroCrossingRateDiff=function(fa, ctrl){
  ZCR=ZerroCrossingRate(fa, ctrl)
  ZCR_1diff=abs(diff(c(0,ZCR)))
  attr(ZerroCrossingRateDiff,"name")="ZCR.D"
  return(ZCR_1diff)
}

ZerroCrossingRateDiff=odf_function("ZCRate",ZerroCrossingRateDiff)


# Amplitude maximum difference
AmplitMaxDiff=function(fa,ctrl) {
  Max.of.wind=apply(fa$ampl, 2, function(x) max(abs(x)))
  AMD=diff(c(0, Max.of.wind))
  return(AMD)
}
AmplitMaxDiff=odf_function("AMax",AmplitMaxDiff)


# Amplitude maximum absolute difference
AmplitMaxAbsDiff=function(fa,ctrl) {
  Max.of.wind=apply(fa$ampl, 2, function(x) max(abs(x)))
  AMD=abs(diff(c(0, Max.of.wind)))
  return(AMD)
}
AmplitMaxAbsDiff=odf_function("AMaxA",AmplitMaxAbsDiff)


# Amplitude energy
AmplitEnergy=function(fa,ctrl) {
  AE=apply(fa$ampl, 2, function(x) sum(x^2))
  return(AE)
}

# Amplitude energy difference
AmplitEnergyDiff=function(fa,ctrl) {
  AE=AmplitEnergy(fa,ctrl)
  AED=diff(c(0, AE)) 
  return(AED)
}
AmplitEnergyDiff=odf_function("AEnerg",AmplitEnergyDiff)


# Amplitude energy absolute difference
AmplitEnergyAbsDiff=function(fa,ctrl) {
  AE=AmplitEnergy(fa,ctrl)
  AED=diff(c(0, AE)) 
  return(abs(AED))
}
AmplitEnergyAbsDiff=odf_function("AEnergA",AmplitEnergyAbsDiff)




# Spectral energy
SpecEnergy=function(fa, ctrl) {
  S = fa$spec
  N=dim(S)[1]
  #S=cbind(rep(0,N), S, rep(0, N)) 
  SE=apply(S, 2, function(x) sum(abs(x)^2))/N
  return(SE)
}

# Spectral energy difference
SpecEnergyDiff=function(fa, ctrl) {
  SE=SpecEnergy(fa, ctrl)
  SE.D=diff(c(0, SE)) 

  return(SE.D)
}
SpecEnergyDiff=odf_function("SEnerg",SpecEnergyDiff)


# Spectral energy absolute difference
SpecEnergyAbsDiff=function(fa, ctrl) {
  SE.D=SpecEnergyDiff(fa, ctrl)
  return(abs(SE.D))
}
SpecEnergyAbsDiff=odf_function("SEnergA",SpecEnergyAbsDiff)



# High frequency content
HighFreqCont=function(fa, ctrl) {
  S = fa$spec
  N=dim(S)[1]
 # S=cbind(rep(0,N), S, rep(0, N)) 
  k=1:N
  HFC=apply(S, 2, function(x) sum(k*abs(x)^2))/N
  return(HFC)
}

# High frequency content difference
HighFreqContDiff=function(fa, ctrl) {
  HFC=HighFreqCont(fa, ctrl)
  HFC.D=diff(c(0, HFC)) 
  return(HFC.D)
}
HighFreqContDiff=odf_function("HFCont",HighFreqContDiff)


# High frequency content absolute difference
HighFreqContAbsDiff=function(fa, ctrl) {
  HFC.D=HighFreqContDiff(fa, ctrl) 
  return(abs(HFC.D))
}
HighFreqContAbsDiff=odf_function("HFContA",HighFreqContAbsDiff)


# Gauss frequency content
HighFreqContGaus=function(fa, ctrl) {
  S = fa$spec
  N=dim(S)[1]
  k=gausswin(N)
  HFCG=apply(S, 2, function(x) sum(k*abs(x)^2))/N
  return(HFCG)
}


# Gauss frequency content difference
HighFreqContGausDiff=function(fa, ctrl) {
  HFCG.D=HighFreqContGaus(fa, ctrl)
  return(diff(c(0,HFCG.D)))
}
HighFreqContGausDiff=odf_function("GFCont",HighFreqContGausDiff)


# Gauss frequency content absolute difference
HighFreqContGausAbsDiff=function(fa, ctrl) {
  HFCG.D=HighFreqContGausDiff(fa, ctrl)
  return(abs(HFCG.D))
}
HighFreqContGausAbsDiff=odf_function("GFContA",HighFreqContGausAbsDiff)



# Spectral centroid
SpecCentroid=function(fa, ctrl){
  S = fa$spec
  k= dim(S)[1] # nr. of frequency bins
  upp = apply(S,2,function(x) {sum(abs(x)*(1:k))})
  low = apply(S,2,function(x) {sum(abs(x))})
  SC = upp / low
  return(SC)
}

# Spectral centroid differenz
SpecCentroidDiff=function(fa, ctrl){
  SC=SpecCentroid(fa, ctrl)
  SC_1diff=abs(diff(c(0,SC)))
  return(SC_1diff)
}
SpecCentroidDiff=odf_function("SCentr",SpecCentroidDiff)


# Spectral spread
SpecSpread=function(fa, ctrl){
  S = fa$spec
  k = dim(S)[1] # nr. of frequency bins
  m = dim(S)[2] # nr. of windows
  SC=SpecCentroid(fa, ctrl)
  
  SS=numeric(m)
  for (i in 1:m) {
    abs.S=abs(S[,i])
    SS[i]=sqrt( sum(abs.S*((1:k)-SC[i])^2)) /
      sqrt(sum(abs.S))
  }
return(SS)
}


# Spectral spread differenz
SpecSpreadDiff=function(fa, ctrl){
  SSp=SpecSpread(fa, ctrl)
  SSp_1diff=abs(diff(c(0,SSp)))
  return(SSp_1diff)
}
SpecSpreadDiff=odf_function("SSpread",SpecSpreadDiff)


# Spectral skewness
SpecSkewness=function(fa, ctrl){
  S = fa$spec
  k = dim(S)[1] # nr. of frequency bins
  m = dim(S)[2] # nr. of windows
  SC=SpecCentroid(fa, ctrl)
  SSp=SpecSpread(fa, ctrl)
  
  SSk=numeric(m)
  for (i in 1:m) {
    abs.S=abs(S[,i])
    SSk[i]=sum(abs.S*((1:k)-SC[i])^3) /
      (SSp[i]^3*(sum(abs.S)))
  }
  return(SSk)
}

# Spectral skewnes difference
SpecSkewnessDiff=function(fa, ctrl){
  SSk=SpecSkewness(fa, ctrl)
  SSk_1diff=abs(diff(c(0,SSk)))
  return(SSk_1diff)
}

SpecSkewnessDiff=odf_function("SSkew",SpecSkewnessDiff)



# Spectral Flux
SpecFlux=function(fa, ctrls) 
{
  S=fa$spec
  n=dim(S)[2]
  H=function(x){
    (x+abs(x))/2
	}
  SF=numeric(n)
  SF[1]=sum(H(abs(S[,1])))
  for (i in 2:n){
    SF[i]=sum( H(abs(S[,i])-abs(S[,(i-1)])) ) 
  }
  return(SF) # damit feauture vector genauso viele Eintr?ge hat wie #Fenster
}
SpecFlux=odf_function("SFlux",SpecFlux)


# Spectral euclid distance
SpecFluxEuclid=function(fa, ctrls)
{
  S=fa$spec
  n=dim(S)[2]
  
  H=function(x){
    x^2
  }
  SF=numeric(n)
  SF[1]=sum(H(abs(S[,1])))
  for (i in 2:n)
  {
    SF[i]=sum( H(abs(S[,i])-abs(S[,(i-1)])) ) 
  }
  return(SF)
}
SpecFluxEuclid=odf_function("SFluxEucl",SpecFluxEuclid)


# Phase deviation
PhaseDeviation=function(fa, ctrl) 
{
  S = fa$spec
  n=dim(S)[2]	# number of windows
  N=dim(S)[1] # number of samples
  S=cbind(rep(0,N),rep(0,N),S)
  
  # intialization
  phi=Arg(S)
  phi_2=phi
  PD=numeric(n+2)
  
  for (i in 3:(n+2)) {
    #phi_1[,i]=phi[,i]-phi[,(i-1)]
    phi_2[,i]=phi[,i]-2*phi[,(i-1)]+phi[,(i-2)]
    PD[i]=(1/dim(S)[1])*sum(abs(phi_2[,i]))
  }
  return(PD[-(1:2)])
}
PhaseDeviation=odf_function("PDev",PhaseDeviation)


# Normalized Weighted phase deviation
WeightPhaseDeviation=function(fa, ctrl) {
  S = fa$spec
  n=dim(S)[2]  # number of windows
  N=dim(S)[1] # number of samples
  S=cbind(rep(0,N),rep(0,N),S)
  
  # intialization
  phi=Arg(S)
  phi_2=phi
  WPD=numeric(n+2)
  
  
  for (i in 3:(n+2)) # for all windows
  {
    #phi_1[,i]=phi[,i]-phi[,(i-1)]
    phi_2[,i]=phi[,i]-2*phi[,(i-1)]+phi[,(i-2)]
    WPD[i]=(1/dim(S)[1])*sum(abs(S[,i]*phi_2[,i]))
  }
  return(WPD[-(1:2)])
}
WeightPhaseDeviation=odf_function("WPDev",WeightPhaseDeviation)


# Complex Domain
ComplexDomain=function(fa, ctrl) {
  S = fa$spec
  n=dim(S)[2]  # number of windows
  N=dim(S)[1] # number of samples

  S=cbind(rep(0,N),rep(0,N),S)
  S.estim=S
  # intialization
  phi=Arg(S)
  phi.est=phi
  CD=numeric(n+2)
  
  for (i in 3:(n+2)) {
    phi.est[,i]=2*phi[,(i-1)]-phi[,(i-2)]
    S.estim[,i]=complex(modulus = Mod(S[,(i-1)]), argument = phi.est[,i])
    CD[i]=(1/dim(S)[1])*sum(abs(S[,i]-S.estim[,i]))
  }
  return(CD[-(1:2)])
}
ComplexDomain=odf_function("CDom",ComplexDomain)


# Rectified Complex Domain
RectComplexDomain=function(fa, ctrl) {
  S = fa$spec
  n=dim(S)[2]  # number of windows
  N=dim(S)[1] # number of samples
  
  S=cbind(rep(0,N),rep(0,N),S)
  S.estim=S
  # intialization
  phi=Arg(S)
  phi.est=phi
  RCD=numeric(n+2)
  
  
    for (i in 3:(n+2)){
    phi.est[,i]=2*phi[,(i-1)]-phi[,(i-2)]
    S.estim[,i]=complex(modulus = Mod(S[,(i-1)]), argument = phi.est[,i])
  }
  
  a.S = abs(S)-abs(S.estim)
  w.pos = a.S > 0
  A = abs(S-S.estim) * w.pos
  RCD = apply(A, 2, sum) * (1/dim(S)[1])
  return(RCD[-(1:2)])
}

RectComplexDomain=odf_function("RCDom",RectComplexDomain)


getODF=function(ctrl){
  switch(ctrl$odf,
         ZCRate=ZerroCrossingRateDiff,
         AMax=AmplitMaxDiff,
         AMaxA=AmplitMaxAbsDiff,
         AEnerg=AmplitEnergyDiff,
         AEnergA=AmplitEnergyAbsDiff,
         SEenerg=SpecEnergyDiff,
         SEenergA=SpecEnergyAbsDiff,
         HFCont=HighFreqContDiff,
         HFContA=HighFreqContAbsDiff,
         GFCont=HighFreqContGausDiff,
         GFContA=HighFreqContGausAbsDiff,
         SCentr=SpecCentroidDiff,
         SSpread=SpecSpreadDiff,
         SSkew=SpecSkewnessDiff,
         SFlux=SpecFlux,
         SFluxEucl=SpecFluxEuclid,
         PDev=PhaseDeviation,
         WPDev=WeightPhaseDeviation,
         CDom=ComplexDomain,
         RCDom=RectComplexDomain,
         # default: try to match the fun which is given as string
         match.fun(ctrl$odf)
  )
}


