[R] fft inverse display help

Dieter Menne dieter.menne at menne-biomed.de
Wed Oct 1 09:24:05 CEST 2008


 <rkevinburton <at> charter.net> writes:

> 
> I have a a simple function that generates a time series square wave:
> 

Your example is not self-running, because the definition of e is unclear.

> Now I ge the x-asis as the real component and y-axis as imaginary component.
When does the display switch?
> How many frequencies do I select before I get the 'Time' axis back? From the
above I am selecting 40 (out of 1500).
> 

When you modify components in the frequency domain, you must always take into
account that these are mirrored at the center coefficient; any changes must be
done symmetrically to the center. Well, strictly speaking: these are periodic,
and it is much easier to think of swapping the upper half in mind to the
negative axis. And a caveat: index 1 (R habits are bad here, it is effectively
frequency 0) is special, it turns up only once and is the DC-component in
electric speak.

Asssume you have 8 coefficients, and you want to lowpass these:

9 4 5 6 1 6 5 4
9 4 5 0 0 0 5 4

As a check if you have done right, look at the imaginary components after the
transformation back to time domain: these should be zero up to rounding errors.

Dieter

genseq <- function()
{
  x <- numeric(4*365)
  s <- seq(as.Date("2005-01-01"), as.Date("2008-12-31"), by="month")
  ob <- as.vector(s[c(10,22,34,46)] - as.Date("2005-01-01"))
  oe <- as.vector(s[c(11,23,35,47)] - as.Date("2005-01-01"))
  for(.index in 1:length(ob))
  {
    x[ob[.index]:(oe[.index] - 1)] <- 1 
  }
  return(ts(x, frequency=365))
}

par(mfcol=c(3,1))
series <- genseq()
fs <- fft(series)
fi <- fft(fs, inverse=TRUE) / length(fs)
plot(fi)

ff <- fs
lowpass <- 10
ff[lowpass:(length(ff)-lowpass+2)] = 0
plot(ff)
fi <- fft(ff, inverse=TRUE) / length(fs)
var(Im(fi)) # Check imaginary part, simple way
plot(fi)

# This contains an error
ff[lowpass:length(ff)] = 0
fi <- fft(ff, inverse=TRUE) / length(fs)
var(Im(fi)) # Imaginary part is large



More information about the R-help mailing list