如何去除时间序列中的特定谐波？

Question

我有一个时间序列，有2014年到2019年的每小时（移动平均）数据。

我使用了 fft() 函数来寻找谐波，以及 quantmod::findPeaks() 找出最高值的位置，结果我得到了以下位置。8 56 2193 4384 6575 8766 306737 308928 311119 313310 315447 315495 (时间序列中共有315499个元素)

现在我需要将其转换回时域，但在这些位置上用零代替数据值。我如何才能做到这一点？

另外，我注意到，当我绘制的是 fft() 函数（FFT的绝对值），x轴还是显示2014年到2019年。是不是应该变成Hz轴（频率）？

z是我的时间序列

FFT <- fft(z)
FFT <- abs(FFT/FFT[1])
plot(FFT)

PK <- findPeaks(FFT, thresh = T)
PK

的结果 PK 是我提到的12个峰值，代表最高的谐波失真.我需要去除这些。我相信我可以简单地做到 FFT[8]=0 到12峰，但我不确定。(我试过，但它不会让我倒退到时域)

我还需要证明，得到的FFT图确实是傅里叶频谱，但x轴不是频率。有什么办法可以进行这种修正吗？

Answer 1

模拟方波+（低）噪声

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并找到峰值

f <- fft(z)
F <- abs(f/f[1])
PK <- quantmod::findPeaks(F,thresh=0.01)

绘制sqrt(功率谱)。

plot(F,type="l")
abline(v=PK,col=2)

剔除峰值

fz <- f
## position of 'peaks' is actually *one before* PK (not sure why?)
fz[PK-1] <- 0

反向变换，与原图比较

zz <- Re(fft(fz,inverse=TRUE))/length(f)
plot(z,type="l")
lines(zz,type="l",col=2)