reconstruction

DESCRIPTION

 Package: gtseries
 Type: Package
-Title: Time Series analysis for paleoclimate and cyclostrat purposes
+Title: Time Series Analysis for Paleoclimate and Cyclostratigraphy Purposes
 Version: 1.05
-Date: 2023-26-10
+Date: 2023-12-15
 Author: Michel Crucifix [aut, cre] (<https://orcid.org/0000-0002-3437-4911>)
 Maintainer: Michel Crucifix <michel.crucifix@uclouvain.be>
 Description: Collection of time series analysis tools, focused on spectral decomposition and reconstruction, with an original focus on astronomical theory of paleoclimates and cyclostratigraphy. 
 License: MIT licence
 LazyLoad: yes
 RoxygenNote: 7.2.3
-Imports: Rdpack, sfsmisc
+Imports: Rdpack, sfsmisc, cmna, gsignal, pracma
 RdMacros: Rdpack
+Depends: 
+    R (>= 2.10)
+LazyData: true

NAMESPACE

@@ -2,6 +2,7 @@
 
 S3method(plot,SSAObject)
 S3method(plot,memObject)
+S3method(plot,mfft_deco)
 S3method(plot,wavelet)
 export(approx_ts)
 export(arspec)
@@ -13,6 +14,7 @@ export(mfft_real)
 export(mfft_real_C)
 export(mfft_real_ter)
 export(powerspectrum.wavelet)
+export(reconstruct_mfft)
 export(ssa)
 importFrom(Rdpack,reprompt)
 importFrom(cmna,goldsectmax)
@@ -23,7 +25,7 @@ importFrom(graphics,mtext)
 importFrom(graphics,par)
 importFrom(graphics,text)
 importFrom(gsignal,hann)
-importFrom(pracma,GramSchmidt)
+importFrom(pracma,gramSchmidt)
 importFrom(sfsmisc,axTexpr)
 importFrom(sfsmisc,eaxis)
 importFrom(stats,ar)

R/mfft.R

@@ -32,8 +32,7 @@
 #' \insertRef{sidlichovsky97aa}{gtseries}
 #
 #' @export mfft
-mfft <- function(xdata, minfreq=NULL, maxfreq=NULL, flag=1, nfreq=30)
-{
+mfft <- function(xdata, minfreq=NULL, maxfreq=NULL, flag=1, nfreq=30) {
 
   xdata = stats::as.ts(xdata)
   dt = deltat(xdata)
@@ -87,10 +86,10 @@ mfft <- function(xdata, minfreq=NULL, maxfreq=NULL, flag=1, nfreq=30)
 
      # if this is a real vector, there will be positive and negative frequencies corresponding to the same amplitude
      # take care and verify that this actually works this way
-
+     OUT = data.frame(Freq=Freq, Ampl=Ampl, Phase=Phase)
+     attr(OUT,"class") = "mfft_deco"
      return(data.frame(Freq=Freq, Ampl=Ampl, Phase=Phase))
-
-  }
+}
 #
 
 #

R/mfft_real.R

@@ -54,8 +54,10 @@ calc_amp <- function(x, nfreq, Fs){
   Freqs <- c(0,Fs)
   amps <- c(sol[1], sqrt ( apply(sol1^2, 2, sum)  ))
   Phases <- -c(0, ( apply(sol1, 2, function(i) Arg(i[1] + 1i*i[2])  )))
-  
-  OUT = data.frame(freqs=Freqs, amps=amps, Phases=Phases) 
+ 
+
+  OUT = data.frame(Freq=Freqs, Amp=amps, Phase=Phases) 
+  attr(OUT,"class") = "mfft_deco"
   return(OUT)
 }
 
@@ -72,7 +74,7 @@ calc_amp <- function(x, nfreq, Fs){
 #' A C-version should be supplied one day. 
 #'
 #' @importFrom cmna goldsectmax
-#' @importFrom pracma GramSchmidt
+#' @importFrom pracma gramSchmidt
 #' @importFrom gsignal hann
 #' @param xdata The data provided either as a time series (advised), or as a vector. 
 #' @param nfreq is the number of frequencies returned, must be smaller that the length of  xdata.
@@ -93,6 +95,7 @@ mfft_real <- function(xdata, nfreq=5){
   dt = deltat(xdata)
   startx = stats::start(xdata)[1]
   N = length(xdata)
+  times <- (seq(N)-1)*dt+startx
   residu = xdata
 
   # will withold the definitive frequencies
@@ -100,14 +103,14 @@ mfft_real <- function(xdata, nfreq=5){
     k=1
     OUT <- calc_amp(xdata, nfreq, Fs)
 
-    print('OUT before corr')
-    print(OUT)
-    Fs <- OUT$freqs[-1]
-    Freqs <- OUT$freqs
-    amps <- OUT$amps
-    Phases <- OUT$Phases
-    # correction 
+    Fs <- OUT$Freq[-1]
 
+
+##   correction  METHOD ANALYTICAL. CURRENTLY NOT USED BUT SHOULD STILL BE EXPLORED
+##  -------------------------------------------------------------------------------- 
+##    Freqs <- OUT$Freq
+##    amps <- OUT$Amp
+##    Phases <- OUT$Phase
 ##     if (k < (nfreq)){
 ##       tau = (N-1) / 2.
 ##       Qp <- function(y) {
@@ -136,18 +139,75 @@ mfft_real <- function(xdata, nfreq=5){
 #   }
    xdata_synthetic <- rep(0,N)
      Fs <- rep(-1, nfreq)
-    for (i in seq(nfreq+1)) xdata_synthetic = xdata_synthetic + OUT$amps[i]*cos(OUT$freqs[i]*seq(N) + OUT$Phases[i])
+    for (i in seq(nfreq+1)) xdata_synthetic = xdata_synthetic + OUT$Amp[i]*cos(OUT$Freq[i]*seq(N) + OUT$Phase[i])
     OUT2 <- calc_amp(xdata_synthetic, nfreq, Fs)
   
-    print(OUT2)
-    print(OUT$freqs - OUT2$freqs)
+#     print(OUT2)
+#     print(OUT$Freq - OUT2$Freq)
     # adjust to units
-    OUT$freqs = OUT$freqs + (OUT$freqs - OUT2$freqs)
-    OUT$amps = OUT$amps + (OUT$amps - OUT2$amps)
-    OUT$Phases = OUT$Phases + (OUT$Phases - OUT2$Phases)
-    OUT$freqs <- OUT$freqs/dt
-    OUT$Phases <- OUT$Phases - startx*OUT$freqs
+    OUT$Freq = OUT$Freq + (OUT$Freq - OUT2$Freq)
+    OUT$Amp = OUT$Amp + (OUT$Amp - OUT2$Amp)
+    OUT$Phase = OUT$Phase + (OUT$Phase - OUT2$Phase)
+    OUT$Freq <- OUT$Freq/dt
+    OUT$Phase <- OUT$Phase - startx*OUT$Freq
+
+    O <- order(OUT$Amp, decreasing=TRUE)
+
+    OUT$Amp = OUT$Amp[O]
+    OUT$Freq = OUT$Freq[O]
+    OUT$Phase = OUT$Phase[O]
+
+  
+
+    attr(OUT,"class") = "mfft_deco"
+    attr(OUT,"nfreq") = nfreq
+    attr(OUT,"xdata") = xdata
+
+    total_ssq <- sum(xdata^2)
+
     return(OUT)
 }
 
+#' MFFT reconstruction
+#' @rdname mfft_deco
+#' @export reconstruct_mfft
+reconstruct_mfft  <- function(M){
+ if (!(attr(M,"class")=="mfft_deco")) stop ("object is not a MFFT decomposition")
+ xdata <- attr(M,"xdata")
+ nfreq <- attr(M,"nfreq")
+ times <- seq(length(xdata))*dt(xdata) + startx(xdata)
+ reconstructed <- lapply(seq(nfreq), function(i) M$Amp[i]*cos(M$Freq[i]*times + M$Phase[i]) )
+}
+
+#' MFFT ANOVA
+#' @rdname mfft_deco
+#' @export mfft_anova
+reconstruct_mfft  <- 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
+plot.mfft_deco <- function (M,...){
+  O <- order(M$Freq)
+  plot(M$Amp[O], M$Freq[O],'h',...)
+}
+
+
 

src/fmft_real.c → deprecated/fmft_real.c

@@ -47,13 +47,13 @@ X(t) + iY(t) = Sum_j=1^N [ A_j * exp i (f_j * t + psi_j) ] */
 #define TWOPI (2.*PI)
 
 
-static float ftemp;
-static double dtemp;
 
 #define DSQR(a) ((dtemp=(a)) == 0.0 ? 0.0 : dtemp*dtemp)
 #define SHFT3(a,b,c) (a)=(b);(b)=(c)
 #define SHFT4(a,b,c,d) (a)=(b);(b)=(c);(c)=(d)
 
+static float ftemp; 
+static double dtemp;
 bool fastflag;
 
 /*bool isPowerofTwo (size_t n){*/