fftを使用したバンドパスフィルタR

Question

私は、サンプリング周波数がfs = 12（月データ）の時系列zを持っており、10月15ヶ月でfftを使ってバンドパスフィルタを実行したいと思います。これは私が進行する方法です：

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で、2つのカラムを持っています。furは高速フーリエ変換（Mod(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はと高くなければなりません.freq - 周波数間隔の境界。詳細は、パッケージドキュメントおよびoctave reference guideに記載されています。

また、fftでは（サンプル周波数）/ 2までの周波数しか得られないことに注意してください。