Perform fft of Time Series in R

Question

I would like to fit waves to a time series using FFT. The goal is to have plots with different harmonics and also use it to forecast n numbers of datapoints.

The code that I'm using is based on this answer from @catastrophic-failure

 nff = function(y = NULL, n = NULL, up = 10L, plot = TRUE, add = FALSE, main = NULL, ...){
  #The direct transformation
  #The first frequency is DC, the rest are duplicated
  dff = fft(y)
  #The time
  t = seq(from = 1, to = length(y))
  #Upsampled time
  nt = seq(from = 1, to = length(y)+1-1/up, by = 1/up)
  #New spectrum
  ndff = array(data = 0, dim = c(length(nt), 1L))
  ndff[1] = dff[1] #Always, it's the DC component
  if(n != 0){
    ndff[2:(n+1)] = dff[2:(n+1)] #The positive frequencies always come first
    #The negative ones are trickier
    ndff[length(ndff):(length(ndff) - n + 1)] = dff[length(y):(length(y) - n + 1)]
  }
  #The inverses
  indff = fft(ndff/as.integer(length(y)), inverse = TRUE)
  idff = fft(dff/as.integer(length(y)), inverse = TRUE)
  if(plot){
    if(!add){
      plot(x = t, y = y, xlab = "Time", ylab = "Data",
           main = ifelse(is.null(main), paste(n, "harmonics"), main), type="l", col="green")
      lines(y = Mod(idff), x = t, col = "red")
    }
    lines(y = Mod(indff), x = nt, ...)
  }
  ret = data.frame(time = nt, y = Mod(indff))
  return(ret)
}

The problem for me is, since I have also negative values in my dataset, that I can´t figure out why positive values are included.

This is the plot of the original data

Compared to the plot after the fft

How can I adapt the code such that the harmonics also cover the missing negative values and also how to use that to calculate (forecast) the next n time points?

Answer 1

The problem appears as you are trying to plot the result with Mod(idff) and Mod(indff) on the following lines:

...
  lines(y = Mod(idff), x = t, col = "red")
}
lines(y = Mod(indff), x = nt, ...)

Mod will always return a positive number corresponding to the magnitude of a complex number.

Since you are computing the inverse FFT on a sequence with Hermitian symmetry (by construction), you should expect a real-valued result. In practice however there may be a small imaginary part due to round-off errors. These may be safely ignored by extracting only the real-parts with Re(idff) and Re(indff), as follows:

...
  lines(y = Re(idff), x = t, col = "red")
}
lines(y = Re(indff), x = nt, ...)

Note that it is usually a good practice to first confirm that the imaginary parts indeed amounting to very small numbers compared with the real parts, since the opposite would suggest that the frequency-domain valued do not have the expected Hermitian symmetry.