user3910073 user3910073 - 8 months ago 59
R Question

Bandpassfilter R using fft

I have a time series

with sampling frequeny
fs = 12
(monthly data) and I would like to perform a bandpass filter using the
at 10 months and 15 months. This is how I would proceed:

y <-
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...

However, I don't know how to derive the frequencies of the fft and I don't know how to plot zz in the time domain. Can someone help me?

Answer Source

I have a function, that wraps fft() a bit:

    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 is values of your signal, and samp.freq is it's sample frequency. It's output is data.frame with two columns - fur is complex numbers we get after fast fourier transform (Mod(fur) will be an amplitude, Arg(fur) - a phase) and freq is vector of corresponding frequencies.

But for frequency filtering I highly reccomend using signal package.

For example using Butterworth filter:

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

In this case interval should be defined as c(Low.freq, High.freq) * (2/samp.freq), where Low.freq and High.freq - borders of frequency intervals. More information can be found in package documentation and octave reference guide.

Also, notice that with fft you can get only frequencies up to (sample frequency)/2.