R help - Jan 2007 - Simple spectral analysis

Hello world,

I am actually trying to transfer a lecture from Statistica to
R and I ran into problems with spectral analysis, I think I
just don't get it 8-(
(The posting from "FFT, frequs, magnitudes, phases" from 2005
did not enlighten me)

As a starter for the students I have a 10year data set of air 
temperature with daily values  and I try to
get a periodogram where the annual period (365 days) should be clearly
visible (in statistica I can get the frequencies and the period).
I tried the spectrum() and pgram() functions, but
did not find a way through... The final aim would be to
get the periodogram (and the residuals from the reassembled data set...)

Thanks and greetings,
Georg

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)


-- 
Georg Hoermann, Dep. of Hydrology, Ecology, Kiel University, Germany
+49/431/23761412, mo: +49/171/4995884, icq:348340729, skype: ghoermann

"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

Hi

without beeing specific in spectrum analysis you will get frequencies 
and spectral densities fro spectrum()
>From help page
An object of class "spec", which is a list containing at least the 
following components: 

freq vector of frequencies at which the spectral density is 
estimated. (Possibly approximate Fourier frequencies.) The units are 
the reciprocal of cycles per unit time (and not per observation 
spacing): see Details below. 
spec Vector (for univariate series) or matrix (for multivariate 
series) of estimates of the spectral density at frequencies 
corresponding to freq. 

<snip>

This is the important part:

**The result is returned invisibly if plot is true.**

So if you call

spectrum(data) you will get plot but in case

sp <- spectrum(data)

you will get also object sp which has above mentioned components. 
Actual periods are obtainable by

n/sp$freq

HTH
Petr

On 8 Jan 2007 at 17:12, Georg Hoermann wrote:

Date sent:      	Mon, 08 Jan 2007 17:12:34 +0100
From:           	Georg Hoermann <georg.hoermann at gmx.de>
To:             	r-help at stat.math.ethz.ch
Subject:        	[R] Simple spectral analysis
> Hello world,
> 
> I am actually trying to transfer a lecture from Statistica to
> R and I ran into problems with spectral analysis, I think I
> just don't get it 8-(
> (The posting from "FFT, frequs, magnitudes, phases" from 2005
> did not enlighten me)
> 
> As a starter for the students I have a 10year data set of air 
> temperature with daily values  and I try to
> get a periodogram where the annual period (365 days) should be clearly
> visible (in statistica I can get the frequencies and the period). I
> tried the spectrum() and pgram() functions, but did not find a way
> through... The final aim would be to get the periodogram (and the
> residuals from the reassembled data set...)
> 
> Thanks and greetings,
> Georg
> 
> 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)
> 
> 
> -- 
> Georg Hoermann, Dep. of Hydrology, Ecology, Kiel University, Germany
> +49/431/23761412, mo: +49/171/4995884, icq:348340729, skype: ghoermann
> 
> ______________________________________________
> R-help at 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.
Petr Pikal
petr.pikal at precheza.cz