[R] Simple spectral analysis

[Earl F. Glynn]

"Georg Hoermann" <georg.hoermann at gmx.de> wrote in message 
news:45A26D72.1040404 at gmx.de...
> The data set:
>
> air = read.csv("http://www.hydrology.uni-kiel.de/~schorsch/air_temp.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("http://www.hydrology.uni-kiel.de/~schorsch/air_temp.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),
     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