r - 使用 fft 的带通滤波器 R

Question

我有一个时间序列z采样频率fs = 12 (每月数据)，我想使用 fft 执行带通滤波器10个月和15个月时。这就是我将如何进行:

y <- as.data.frame(fft(z))
y$freq <- ..
y$y <-  ifelse(y$freq>= 1/10 & y$freq<= 1/15,y$y,0)
zz <- fft(y$y, inverse = TRUE)/length(z)
plot zz in the time domain...

但是，我不知道如何导出 fft 的频率，也不知道如何在时域中绘制 zz。有人可以帮助我吗？

Answer 1

我有一个函数，它包装了 fft() 一点:

    function(y, samp.freq, ...){
      N <- length(y)
      fk <- fft(y)
      fk <- fk[2:length(fk)/2+1]
      fk <- 2*fk[seq(1, length(fk), by = 2)]/N
      freq <- (1:(length(fk)))* samp.freq/(2*length(fk))
      return(data.frame(fur = fk, freq = freq))
    }

y 是信号值，samp.freq 是它的采样频率。它的输出是带有两列的 data.frame - fur 是我们在快速傅里叶变换后得到的复数(Mod(fur) 将是一个幅度，Arg(fur) - 相位)，freq 是相应频率的向量。

但是对于频率过滤，我强烈建议使用信号包。

例如使用巴特沃斯滤波器:

     library('signal')
     bf <- butter(2, c(low, high), type = "pass")
     signal.filtered <- filtfilt(bf, signal.noisy)

在这种情况下，间隔应定义为 c(Low.freq, High.freq) * (2/samp.freq)，其中 Low.freq 和 High.freq - 频率间隔的边界。更多信息可以在包文档和 octave reference guide 中找到。 .

此外，请注意，使用 fft 只能获得高达(采样频率)/2 的频率。