Newer
Older
#' @param M : mfft_deco object
#' @param sum : TRUE if user wants to sum components in the reconstruction
#' @return list of reconstructed components if sum=FALSE, full
#' reconstructed time series otherwise
reconstruct_mfft <- function(M, sum=TRUE){
if (!(attr(M,"class")=="mfft_deco")) stop ("object is not a MFFT decomposition")
times <- seq(length(xdata))*stats::deltat(xdata) + stats::start(xdata)[1]
reconstructed <- lapply(seq(nfreq), function(i) ts( M$Amp[i]*cos(M$Freq[i]*times + M$Phase[i]), start=stats::start(xdata), deltat=stats::deltat(xdata)) )
if ( sum ) reconstructed <- ts(apply(simplify2array(reconstructed), 1 , sum), start=stats::start(xdata), deltat=stats::deltat(xdata))
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
}
#' MFFT ANOVA
#' not ready. do not use.
#' @rdname mfft_deco
#' @export mfft_anova
mfft_anova <- function(M){
if (!(attr(M,"class")=="mfft_deco")) stop ("object is not a MFFT decomposition")
xdata <- attr(M,"xdata")
nfreq <- attr(M,"nfreq")
N <- length(xdata)
times <- seq(length(xdata))*deltat(xdata) + start(xdata)
reconstructed <- sapply(seq(nfreq), function(i) M$Amp[i]*cos(M$Freq[i]*times + M$Phase[i]) )
cum_reconstruct <- apply(reconstructed, 1, cumsum)
residual_vars <- apply(apply(cum_reconstruct, 1, function(i) xdata-i) , 2, function(j) sum(j^2))
var0 <- sum(xdata^2)
master_vars <- c(var0, residual_vars[-length(residual_vars)])
p2s <- 2*seq(nfreq)
p1s <- c(0, p2s[-length(p2s)])
F <- (master_vars - residual_vars)/residual_vars * (N - p2s)/(p2s-p1s)
}
#' @rdname mfft_deco
#' @export
as.data.frame.mfft_deco <- function(x) {data.frame(Freq=x$Freq, Amp=x$Amp, Phases=x$Phases)}
#' @rdname mfft_deco
#' @export
plot.mfft_deco <- function (M,periods=FALSE,...){
# O <- order(M$Freq)
plot(abs(M$Freq), abs(M$Amp),'h',ylab="Amplitudes", xlab="", ...)
if (periods) {
frequencies <- pretty(range(M$Freq/(2*pi)))
labels <- as.character(1/frequencies)
if (0 %in% frequencies) labels[which(frequencies == 0)] = "∞"
axis(1, line=3, at=2*pi*frequencies, labels=labels)
mtext("Rate", 1, 2)
mtext("Period", 1, 4)
} else {
mtext("Rate", 1, 3)
}
points(abs(M$Freq), abs(M$Amp),'p',...)
}
#' @rdname mfft_deco
#' @export
lines.mfft_deco <- function (M,...){
# O <- order(M$Freq)
lines(abs(M$Freq), abs(M$Amp),'h',...)
points(abs(M$Freq), abs(M$Amp),'p',...)
}
#' @rdname mfft_deco
#' @export
print.mfft_deco <- function (M,...){
print.data.frame(cbind(as.data.frame(M), Period=2*pi/M$Freq))
}