Re: [R] Simple spectral analysis

From: Earl F. Glynn <>
Date: Mon 08 Jan 2007 - 22:07:29 GMT

"Georg Hoermann" <> wrote in message
> The data set:
> air = read.csv("")
> airtemp = ts(T_air, start=c(1989,1), freq = 365)
> plot(airtemp)

Maybe this will get you started using fft or spectrum -- I'm not sure why the spectrum answer is only close:

air = read.csv("")

TempAirC <- air$T_air
Time <- as.Date(air$Date, "%d.%m.%Y") N <- length(Time)

oldpar <- par(mfrow=c(4,1))
plot(TempAirC ~ Time)

# Using fft
transform <- fft(TempAirC)

# Extract DC component from transform
dc <- Mod(transform[1])/N

periodogram <- round( Mod(transform)^2/N, 3)

# Drop first element, which is the mean
periodogram <- periodogram[-1]

# keep first half up to Nyquist limit
periodogram <- periodogram[1:(N/2)]

# Approximate number of data points in single cycle: print( N / which(max(periodogram) == periodogram) )

# plot spectrum against Fourier Frequency index plot(periodogram, col="red", type="o",

     xlab="Fourier Frequency Index", xlim=c(0,25),
     main="Periodogram derived from 'fft'")

# Using spectrum
s <- spectrum(TempAirC, taper=0, detrend=FALSE, col="red",

              main="Spectral Density")

plot(log(s$spec) ~ s$freq, col="red", type="o",

     xlab="Fourier Frequency", xlim=c(0.0, 0.005),
     main="Periodogram from 'spectrum'")

cat("Max frequency\n")
maxfreq <- s$freq[ which(max(s$spec) == s$spec) ]

# Period will be 1/frequency:
cat("Corresponding period\n")



Earl F. Glynn
Scientific Programmer
Stowers Institute mailing list PLEASE do read the posting guide and provide commented, minimal, self-contained, reproducible code. Received on Tue Jan 09 18:55:57 2007

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Tue 09 Jan 2007 - 08:30:29 GMT.

Mailing list information is available at Please read the posting guide before posting to the list.