# Re: [R] Simple spectral analysis

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

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

R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html 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 https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.