使用FFT计算带通滤波器R

3

我有一个时间序列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。有人能帮我吗?

1个回答

5

我有一个函数,它稍微包装了一下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是相应频率的向量。

但是为了进行频率过滤,我强烈推荐使用signal包。

例如使用Butterworth滤波器:

     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参考指南
另外,请注意,使用fft只能获得高达(采样频率)/2的频率。

谢谢,但这些频率与我的带通窗口(10和15)有什么关系?我如何将不对应于此间隔的fft值设置为零? - user3910073
1
更新了答案,其中包含适用于您的值。 - Andrey Kolyadin

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接