Earl F. Glynn wrote: > "Georg Hoermann" <[EMAIL PROTECTED]> wrote in message > news:[EMAIL PROTECTED] > >> 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.