Alvaro Aballe Villero wrote:
> Let`s consider a vector "ni, i=1,2,..N"...
> If I perfom the a FFT, I expect to obtain another vector "Ck, k=1,2,..N/2"
> ...However, with the function "fft" I get a N-point long result.
The Fourier transform of a *complex* N-vector "z" is another complex N-vector
"w". If "z" is *real*, then the second half of "w" is the complex conjugate of
the first half, so it can be ignored.
For z real, let:
N <- length(z)
w <- fft(z)/N # Division by N here is a little unconventional
jmax <- floor(N/2+1)
tv <- 0:(N-1) # Time values, starting at 0
w[1]==mean(z) is the constant term.
For j=2:jmax, w[j] is the Fourier coefficient for frequency (j-1)/N, and
w[N+2-j] = Conj(w[j]). Note that for even N, this implies w[jmax] is real.
You could derive the Fourier transform yourself (more slowly) with:
w <- rep(NA, N)
w[1] <- mean(z)
for (j in 2:jmax) {
w[j] <- sum(z * exp(-2i*pi*tv*(j-1)/N)) / N
w[N+2-j] <- Conj(w[j])
}
or even:
w <- rep(NA, N)
w[1] <- mean(z)
for (j in 2:jmax) {
phi <- 2*pi*tv*(j-1)/N
w[j] <- sum(z * (cos(phi) - 1i*sin(phi))) / N
w[N+2-j] <- sum(z * (cos(phi) + 1i*sin(phi))) / N
}
You can reconstruct the original series "z" either by:
zz <- Re(fft(w, inverse=T))
or (more slowly) by:
zz <- rep(Re(w[1]), N)
for (j in 2:jmax) {
dup <- if (N%%2==0 && j==jmax) 1 else 2 # Don't double jmax if N even
zz <- zz + dup * Mod(w[j]) * cos(Arg(w[j]) + 2*pi*tv*(j-1)/N)
}
The power spectrum ("periodogram"):
s <- N * abs(w[2:jmax])^2
plot((2:jmax-1)/N, s, type="l", log="y", xlab="frequency", ylab="spectrum")
can also be obtained from the "spectrum" function:
s1 <- spectrum(z, fast=F, taper=0, detrend=F) # Plots automatically
plot(s1$freq, s1$spec, type="l", log="y", xlab="frequency", ylab="spectrum")
Note I had to override the default to detrend, taper, and pad z. Smoother
(more meaningful?) spectra are obtained with, e.g.:
s2 <- spectrum(z, spans=c(2*round(N/20)+1, 2*round(N/10)+1))
s3 <- spec.ar(z)
Among other things, the power spectrum shows the distribution in frequency
space of the "sum-of-squares" errors. Again you must watch that tricky jmax:
sp <- s; if (N%%2==0) sp[jmax-1] <- sp[jmax-1]/2
but then (N-1)*var(z) == 2*sum(sp).
--
-- David Brahm (brahm@alum.mit.edu)
|