分数阶微分 FOD(GL)

分数阶微分方程讲义-华师大数院

分数阶微分-清华大学出版社

参考以上文件公式与代码，使用R代码编写，针对光谱数据，经对比alpha=1时与应用prospectr包函数计算结果一致

fod<- function(spectrum, alpha) {
  # init
  d <- dim(spectrum)
  FOD <- matrix(0, nrow = 0, ncol = d[2])

  # calculate h<band width>
  wavelength12 <- as.integer(unlist(lapply(wavelength12 <- colnames(spectrum)[1:2],function(x) substr(x,2,nchar(x)))))
  h <-wavelength12[2]-wavelength12[1]

  # calculate W
  w <- c(1, rep(0, d[2] - 1))
  for (j in 2:d[2])
    w[j] <- w[j - 1] * (1 - (alpha + 1) / (j - 1))

  # calculate fod
  for (fid in 1:d[1]) {
    thisSpec <- as.numeric(spectrum[fid,])
    thisFod <-  vector(mode='numeric',length=d[2])

    for (i in 1:d[2]) {
      for (j in 0:(i - 1))
        thisFod[i] <- thisFod[i] + thisSpec[i - j] * w[j + 1]
    }
    FOD <- rbind(FOD, thisFod)
    cat(fid,"\t")
  }

  FOD <- FOD / (h ^ alpha)

  # format
  FOD <- as.data.frame(FOD)
  dimnames(FOD) <- dimnames(spectrum)


  #
  FOD = FOD[,-(1:6)]

  cat("FOD ",alpha,"finished\n")
  return (FOD)
}