我有一个z采样频率fs = 12(每月数据)的时间序列,我想使用fft10 个月和 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...
Run Code Online (Sandbox Code Playgroud)
但是,我不知道如何推导出 fft 的频率,也不知道如何在时域中绘制 zz。有人能帮我吗?
我有一个函数,它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))
}
Run Code Online (Sandbox Code Playgroud)
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)
Run Code Online (Sandbox Code Playgroud)
在这种情况下,间隔应定义为 c(Low.freq, High.freq) * (2/samp.freq),其中 Low.freq 和 High.freq - 频率间隔的边界。更多信息可以在包文档和八度参考指南中找到。
另外,请注意,使用 fft 您只能获得高达(采样频率)/2 的频率。
| 归档时间: |
|
| 查看次数: |
1139 次 |
| 最近记录: |