[R] Synthesis of harmonic functions

[Thomas Petzoldt]

Hello,

I try to synthesize harmonic functions with a subset of frequencies
determined by fast fourier transform. I wrote two different functions, a
looped version and a matrix multiplication version. As an example the
looped version takes 5 sec and the matrix-version takes 3-4 sec (R 1.4,
Athlon 1.2 GHz, 256 MB, Win ME), but the latter needs huge amount of
memory due to the matrices. So, none of them is satisfying.

Is there a common way, a better solution or a built-in function (e.g.
inverse of fft) which I have overlooked?

Thank you very much!

Thomas Petzoldt

##====================================================

harmonic.simple <- function(x, a0, a, b, t, ord) {
  y <- a0
  for (p in ord) {
    k <- 2 * pi * p * x/t
    y <- y + a[p] * cos(k) + b[p] * sin(k)
  }
  y
}
   
harmonic.matrix <- function(x, a0, a, b, n, ord) {
  k <- (2 * pi * x/n) %*% t(ord)
  y <- a0 + cos(k) %*% a[ord] + sin(k) %*% b[ord]
  y 
}

n      <- 4000
maxord <- n/2
x      <- 0:(n-1)
y      <- sin(x*2*pi/n) + rnorm(n)*0.1

pf <- fft(y)      ## Fast Fourier Transform
a0 <- Re(pf[1])/n
pf <- pf[-1]      ## drop first element
a  <-  2*Re(pf)/n; b  <- -2*Im(pf)/n

y1 <-  harmonic.simple(x,a0,a,b,n,1:maxord)
y2 <-  harmonic.matrix(x,a0,a,b,n,1:maxord)

#plot(x,y)
#lines(x, y1,col="green")
#lines(x, y2,col="red",lty="dashed")

##====================================================
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._