# [R] Simple spectral analysis

Earl F. Glynn efg at stowers-institute.org
Mon Jan 8 23:07:29 CET 2007

```"Georg Hoermann" <georg.hoermann at gmx.de> wrote in message
news:45A26D72.1040404 at gmx.de...
> The data set:
>
> 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:

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),
ylab="Periodogram",
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),
ylab="Log(Periodogram)",
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")
print(1/maxfreq)

par(oldpar)

efg

Earl F. Glynn
Scientific Programmer
Stowers Institute

```