[R] Simple spectral analysis

[Peter Dalgaard]

Earl F. Glynn wrote:
> "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:
>   
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 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.
>