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