#plotFFT_function
#x: data_time_scale
#y: data
#samplingFreq: the sampling frequence of data capture device
#dominant_frequency_order: defult first (1)
#cex.axis: the size of head note for dominant frequencey
#PSD : Power Spectral Density (modFFT**2[2:(N/2+1)]*2/N)

plotFFT <- function(x, y, samplingFreq, dominant_frequency_order=1,
                    plotFreqs = (samplingFreq/2),
                    cex.axis = 1, PSD = F, shadeNyq=TRUE){
  getFFTFreqs <- function(Nyq.Freq, data){
  if ((length(data) %% 2) == 1) # Odd number of samples
  {
    FFTFreqs <- c(seq(0, Nyq.Freq, length.out=(length(data)+1)/2),
                  seq(-Nyq.Freq, 0, length.out=(length(data)-1)/2))
  } else{
    FFTFreqs <- c(seq(0, Nyq.Freq, length.out=length(data)/2),
                  seq(-Nyq.Freq, 0, length.out=length(data)/2))
  }# Even number

  return (FFTFreqs)
}
  if(PSD == F){
  Nyq.Freq <- samplingFreq/2
  FFTFreqs <- getFFTFreqs(Nyq.Freq, y)
  FFT <- fft(y)
  modFFT <- Mod(FFT)
  FFTdata <- data.frame(cbind(FFTFreqs, modFFT=round(modFFT,3)))
  plot(FFTdata[FFTdata$FFTFreqs>0 & FFTdata$FFTFreqs < plotFreqs ,], t="l", pch=20, lwd=2, cex=0.8, main="",
       xlab="Frequency (Hz)", ylab="Power")

  FFTdata1<-FFTdata[FFTdata$FFTFreqs > 1,]
  dominant_frequency<-round(FFTdata1[order(FFTdata1$modFFT,decreasing=TRUE),]$FFTFreqs[1:5],2)
    # Period axis on top
    axis(3, cex.axis=cex.axis,
         labels = paste0("dominant_frequency:", dominant_frequency[dominant_frequency_order]),
         at=dominant_frequency[dominant_frequency_order])
  if (shadeNyq == TRUE)
  {
    # Gray out lower frequencies
    rect(dominant_frequency[dominant_frequency_order]+0.5, 0, dominant_frequency[dominant_frequency_order]-0.5, max(FFTdata[,2]), col="gray", density=30)
  }
  #ret <- list("freq"=FFTFreqs, "FFT"=FFT, "modFFT"=modFFT)
  return (dominant_frequency = dominant_frequency[dominant_frequency_order])
}
  if(PSD == T){
  #add PSD_112017
  Nyq.Freq <- samplingFreq/2
  Nf = length(y)/2
  FFTFreqs <- seq.int(from=0, to=Nyq.Freq, length.out=Nf)
  FFT <- fft(y)
  modFFT <- Mod(FFT)

  #add PSD_112017
  Se = modFFT**2
  S <- c(Se[2:(Nf+1)] * 2 / Nyq.Freq) # Finally, the PSD!

  #plot PSD
  FFTdata <- data.frame(cbind(FFTFreqs, modFFT=round(S,3)))
  plot(FFTdata[FFTdata$FFTFreqs>0 & FFTdata$FFTFreqs < plotFreqs ,], t="l", pch=20, lwd=2, cex=0.8, main="",
       xlab="Frequency (Hz)", ylab="PSD")
  FFTdata1<-FFTdata[FFTdata$FFTFreqs > 1,]
  dominant_frequency<-round(FFTdata1[order(FFTdata1$modFFT,decreasing=TRUE),]$FFTFreqs[1:5],2)
  # Period axis on top
  axis(3, cex.axis=cex.axis,
       labels = paste0("dominant_frequency:", dominant_frequency[dominant_frequency_order]),
       at=dominant_frequency[dominant_frequency_order])
  if (shadeNyq == TRUE)
  {
    # Gray out lower frequencies
    rect(dominant_frequency[dominant_frequency_order]+0.5, 0, dominant_frequency[dominant_frequency_order]-0.5, max(FFTdata[,2]), col="gray", density=30)
  }
  return (dominant_frequency = dominant_frequency[dominant_frequency_order])}

}