Re: [R] Simple spectral analysis

From: Peter Dalgaard <p.dalgaard_at_biostat.ku.dk>
Date: Mon 08 Jan 2007 - 23:37:58 GMT

Earl F. Glynn wrote:
> "Georg Hoermann" <georg.hoermann@gmx.de> wrote in message
> news:45A26D72.1040404@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:
>
The defaults for detrending and tapering could be involved. Putting, e.g., detrend=F gives me a spectrum with substantially higher low-frequency components.

But what was the problem in the first place?

spec.pgram(airtemp,xlim=c(0,10))

abline(v=1:10,col="red")

shows a strong peak at 1 and maybe a weak peak at 2, and the other integer frequencies less pronounced. This seems reasonably in tune with

> x <- (1:3652)/365

> summary(lm(air$T_air ~ sin(2*pi*x)+cos(2*pi*x)+ sin(4*pi*x)+cos(4*pi*x) + sin(6*pi*x)+cos(6*pi*x)+x))

Call:

lm(formula = air$T_air ~ sin(2 * pi * x) + cos(2 * pi * x) +

    sin(4 * pi * x) + cos(4 * pi * x) + sin(6 * pi * x) + cos(6 *

    pi * x) + x)

Residuals:

     Min 1Q Median 3Q Max

-16.3109 -2.3317 -0.1080 2.2063 10.6249

Coefficients:

                Estimate Std. Error t value Pr(>|t|)    

(Intercept)      9.67904    0.11267  85.909   <2e-16 ***

sin(2 * pi * x) -2.64554 0.07967 -33.208 <2e-16 ***

cos(2 * pi * x) -7.73520 0.07938 -97.443 <2e-16 ***

sin(4 * pi * x) 0.92967 0.07948 11.696 <2e-16 ***

cos(4 * pi * x) 0.13982 0.07938 1.761 0.0783 .

sin(6 * pi * x) 0.13320 0.07945 1.676 0.0937 .

cos(6 * pi * x) 0.14480 0.07938 1.824 0.0682 .

x -0.23773 0.01952 -12.179 <2e-16 ***

---

Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1   1 


Residual standard error: 3.393 on 3644 degrees of freedom

Multiple R-Squared: 0.7486,     Adjusted R-squared: 0.7482 

F-statistic:  1550 on 7 and 3644 DF,  p-value: < 2.2e-16 






> 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
>
> ______________________________________________
> 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.
>
______________________________________________ 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 17:07:46 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 - 10:30:30 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.