I have a time series with hourly (moving average) data from 2014 to 2019.
I used the fft() function to find the harmonics, and quantmod::findPeaks() to find the positions of the highest valued ones, which gave me as result the following positions: 8 56 2193 4384 6575 8766 306737 308928 311119 313310 315447 315495 (there are 315499 elements in the time series).
Now I need to transform this back to the time domain, but substituting zero for the value of the data in those positions. How can I do that?
Also, I noticed that when I plotted the fft() function (the absolute values of the FFT), the x-axis was still showing the years 2014 to 2019. Shouldn't it turn into Hz axis (frequency)?
z is my time series
FFT <- fft(z)
FFT <- abs(FFT/FFT[1])
plot(FFT)
PK <- findPeaks(FFT, thresh = T)
PK
The results of PK are the 12 peaks I mentioned, representing the highest harmonic distortions.
I need to remove those. I believe I can simply do FFT[8]=0 to the 12 peaks, but I'm not sure. (I tried that, but it won't let me reverse back to time domain)
I also need to prove that the resulting FFT plot is really the Fourier spectrum, but the x-axis is not frequency. Any idea on how to make that correction?
Simulated square wave + (low) noise
set.seed(101)
ncycles <- 20
len <- 50
z <- rep(rep(0:1,each=len),ncycles) + rnorm(len*ncycles,mean=0, sd=0.02)
plot(z,type="l")
FFT and find peaks
f <- fft(z)
F <- abs(f/f[1])
PK <- quantmod::findPeaks(F,thresh=0.01)
Plot sqrt(power spectrum):
plot(F,type="l")
abline(v=PK,col=2)
Zero out peaks:
fz <- f
## position of 'peaks' is actually *one before* PK (not sure why?)
fz[PK-1] <- 0
Inverse transform, compare with original:
zz <- Re(fft(fz,inverse=TRUE))/length(f)
plot(z,type="l")
lines(zz,type="l",col=2)