wrap up + doc

NAMESPACE

@@ -13,12 +13,9 @@ export(arspec)
 export(cwt_morlet)
 export(hilbert_extension)
 export(mem)
-export(mfft)
 export(mfft_anova)
+export(mfft_complex)
 export(mfft_real)
-export(mfft_real_C)
-export(mfft_real_quatro)
-export(mfft_real_ter)
 export(periodogram)
 export(powerspectrum.wavelet)
 export(reconstruct_mfft)
@@ -32,8 +29,6 @@ importFrom(graphics,lines)
 importFrom(graphics,mtext)
 importFrom(graphics,par)
 importFrom(graphics,text)
-importFrom(gsignal,hann)
-importFrom(pracma,gramSchmidt)
 importFrom(sfsmisc,axTexpr)
 importFrom(sfsmisc,eaxis)
 importFrom(stats,ar)

R/mfft.R → R/mfft_complex.R

-# I need to add lots of comments about licencing
-# the fact that the hanning window is now made here
-# check that it provides the same output as the original code
-# espacially with the new window
+#' Modified Fourier transform  
+#'
+#' Implementation of the the Frequency Modified Fourier Transform 
+#' (Sidlichovsky and Nesvorny 1997, Cel. Mech. 65, 137). 
+#' Given a quasi-periodic complex signal X + iY, the algorithm 
+#' estimates the frequencies (f_j), amplitudes (A_j) and phases 
+#' (psi_j) in its decomposition. The algorithm was generalised
+#' by M. Crucifix when the input signal is real. 
+#' `mfft` is a wrapper thar redirects the data
+#' to `mfft_real` or `mfft_complex` depending on the nature of
+#' the input data. 
+#' TODO: this needs more work  because the two routines
+#' do not yet take exactly the same inputs. In particular, the meaning
+#' the `correction` parameter is not the same 
+#' TODO: add Laskar reference and say that the correction=0 version
+#'       corresponds to Laskar's routine.
+#' @param xdata The data provided either as a time series (advised), or as a vector. 
+#' may be complex
+#' @param minfreq,maxfreq If provided, bracket the frequencies to be probed. Note this are
+#'        more precisely angular velocities (2\pi / period), expressed in time-inverse units
+#'        with the time resolution encoded in `xdata` if the latter is a time series. 
+#'        The computed frequencies are in the range given by minfreq and maxfreq.
+#' @param correction :
+#'        The parameter is passed on to `mfft_real` or `mfft_complex` depending on
+#'        `xdata` is real or complex, but the meaning of the parameter is not
+#'        quite the same for both. This needs to be fixed. 
+#'        one approach would be to define scalars like 'fft', 'mfft', 'mfft_linear', 
+#'        'mfft_second_order_correction'. 
+#' @param nfreq is the number of frequencies returned, must be smaller that the length of  xdata.
+#' @return a `mfft_deco` object, based on a data.frame with columns "Freq", "Ampl" and "Phases". 
+#'         note that because of a language glitch (to be fixed), "Freq" actually means "Rate"
+#' @author Michel Crucifix for the R code, and David Nesvorny for most of the supporting C code doing the
+#'         Frequency Modified Fourier transform  for Complex Numbers
+#' \insertRef{sidlichovsky97aa}{gtseries}
+mfft <- function(xdata, nfreq=30, minfreq=NULL, maxfreq=NULL, correction=1) {
+  if (is.complex(xdata)) return(mfft_complex(xdata, nfreq, minfreq, maxfreq, correction)) else 
+                         return(mfft_real(xdata, nfreq, minfreq, maxfreq, correction)) 
+}
 
-#' Modified Fourier transform 
+#' Frequency Modified Fourier transform  for Complex Numbers
 #'
 #' Implementation of the the Frequency Modified Fourier Transform 
 #' (Sidlichovsky and Nesvorny 1997, Cel. Mech. 65, 137). 
@@ -18,64 +52,52 @@
 #' @param minfreq,maxfreq If provided, bracket the frequencies to be probed. Note this are
 #'        more precisely angular velocities (2\pi / period), expressed in time-inverse units
 #'        with the time resolution encoded in `xdata` if the latter is a time series. 
-#' @param flag :
-#'      Modified Fourier Transform                  if   flag = 1;
-#'      Frequency Modified Fourier Transform        if   flag = 2;
-#'      FMFT with additional non-linear correction  if   flag = 3
+#' @param correction :
+#'      Modified Fourier Transform                  if   correction = 0;
+#'      Frequency Modified Fourier Transform        if   correction = 1;
+#'      FMFT with additional non-linear correction  if   correction = 2
 #' (while the first algorithm is app. 3 times faster than the third one, 
 #' the third algorithm should be in general much more precise).  
 #' The computed frequencies are in the range given by minfreq and maxfreq.
 #' @param nfreq is the number of frequencies returned, must be smaller that the length of  xdata.
-#' @return STILL NEED TO BE DESCRIBED
+#' @return a `mfft_deco` object, based on a data.frame with columns "Freq", "Ampl" and "Phases". 
+#'         note that because of a language glitch (to be fixed), "Freq" actually means "Rate"
 #' @author Michel Crucifix for the R code, and David Nesvorny for most of the supporting C code doing the
 #' actual computations
 #' \insertRef{sidlichovsky97aa}{gtseries}
 #
-#' @export mfft
-mfft <- function(xdata, minfreq=NULL, maxfreq=NULL, flag=1, nfreq=30) {
+#' @export mfft_complex
+mfft_complex <- function(xdata, nfreq=30,  minfreq=NULL, maxfreq=NULL, correction=1) {
 
-  xdata = stats::as.ts(xdata)
-  dt = deltat(xdata)
-  startx = stats::start(xdata)[1]
-  ydata = Im(xdata)
-  xdata = Re(xdata)
-  ndata = length(xdata);
+  xdata <- stats::as.ts(xdata)
+  dt <- deltat(xdata)
+  startx <- stats::start(xdata)[1]
+  ydata <- Im(xdata)
+  xdata <- Re(xdata)
+  ndata <- length(xdata);
 
+  flag <- correction + 1  # this is an ugly translation towards Sidlichovsky's convention
 
   if (!(flag == 1 || flag == 2 || flag == 3)) stop("flag must be either 1, 2, or 3")
 
-#   if (window=="hanning")
-#   {
-#     hanning = (1-cos(2*pi*(seq(ndata)-1)/(ndata-1)))/2.
-#     windowing
-#     xdata = xdata*hanning
-#     ydata = ydata*hanning
-#   } else if (window=="none") {} else stop ("this window is not supported, sorry")
-#     padding
-#     N = 2^ceiling(log(ndata)/log(2))
-#     if (N > ndata){
-#       xdata[(ndata+1):N] = 0
-#       ydata[(ndata+1):N] = 0
-#     }
-#     
 
     if (is.null(minfreq)){
-      my_minfreq = -pi
-      my_maxfreq = pi} else {
-      my_minfreq = minfreq * dt
-      my_maxfreq = maxfreq * dt}
+     my_minfreq <- -pi
+     my_maxfreq <- pi} else {
+     my_minfreq <- minfreq * dt
+     my_maxfreq <- maxfreq * dt}
  
-     signal1 = as.double(rep(0,3*nfreq))
-     signal2 = as.double(rep(0,3*nfreq))
-     signal3 = as.double(rep(0,3*nfreq))
+     signal1 <- as.double(rep(0,3*nfreq))
+     signal2 <- as.double(rep(0,3*nfreq))
+     signal3 <- as.double(rep(0,3*nfreq))
 
-     Infreq=as.integer(nfreq)
-     Iminfreq = as.double(my_minfreq)
-     Imaxfreq = as.double(my_maxfreq)
-     Ixdata = as.double(xdata)
-     Iydata = as.double(ydata)
+     Infreq<-as.integer(nfreq)
+     Iminfreq <- as.double(my_minfreq)
+     Imaxfreq <- as.double(my_maxfreq)
+     Ixdata <- as.double(xdata)
+     Iydata <- as.double(ydata)
 
-     OUT = .C("fmft", Infreq, Iminfreq, Imaxfreq, as.integer(flag), 
+     OUT <- .C("fmft", Infreq, Iminfreq, Imaxfreq, as.integer(flag), 
               as.integer(ndata), Ixdata, Iydata, signal1,signal2,signal3, DUP=TRUE)
 
      OUTVEC <- t(matrix(OUT[[7+flag]], 3, nfreq))
@@ -84,17 +106,21 @@ mfft <- function(xdata, minfreq=NULL, maxfreq=NULL, flag=1, nfreq=30) {
      Ampl <- OUTVEC[nfreq+seq(nfreq) ] 
      Phase <- OUTVEC[2*nfreq+seq(nfreq) ] - startx*Freq
 
-     # 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, Amp=Ampl, Phases=Phase)
-     attr(OUT,"class") = "mfft_deco"
-     attr(OUT,"nfreq") = nfreq
-     attr(OUT,"xdata") = xdata
+     # if this is input is a  real vector, there will be positive and negative frequencies 
+     # corresponding to the same amplitude
+     # however, better use mfft_real in this case
+     OUT <- data.frame(Freq=Freq, Amp=Ampl, Phases=Phase)
+
+     attr(OUT,"class") <- "mfft_deco"
+     attr(OUT,"nfreq") <- nfreq
+     attr(OUT,"xdata") <- xdata
+
      return(OUT)
 }
 
 
 
+
 #
 
 #

R/mfft_real_quatro.R

-
-Q <- function(wT) {sin(wT)/(wT)*(pi^2)/(pi^2-(wT)^2)}
-Qprime <- function(y) {ifelse(y==0, 0, pi^2/(pi^2-y^2)/y*(cos(y)+(sin(y)/y)*(3*y^2-pi^2)/(pi^2-y^2)))}
-Qsecond0 <- 2/pi^2 - 1./3. 
-
-analyse_quatro <- function(xdata, nfreq, fast = TRUE, nu = NULL){
-  
-
-  # nu = c(0.0429, 0.0186)
-  # nu can be provided, in which case the frequencies
-  # are considered to be given
-  if (!is.null(nu))
-  {
-    # we will need two versions of nfreq, one for the cosine and one for the cosine
-    # if if there is a zero frequency we only need it once
-    nutemp <- unlist ( lapply( nu, function(a) if (a==0) return (a) else return (  c(a,-a)))  )
-    phase <-  unlist ( lapply( nu, function(a) if (a==0) return (0) else return (  c(0, pi/2)))  )
-    nu <- nutemp
-    if (length(nu) < 2*nfreq) {
-      # this will typically occur if one zero frequency was provided
-      nu[2*nfreq] <- NA
-      phase[2*nfreq] <- NA
-    }
-  } else 
-  {
-    nu <- rep(NA,2*nfreq)
-    phase <- rep(NA,2*nfreq)
-  }
-
-  N <- length(xdata)
-  T <- N / 2
-  hann <- function(N) {return (1-cos(2*pi*seq(0,(N-1))/(N)))};
-  hN <- hann(N)
-  t <- seq(N)-1
-  power <- function(x) (Mod(fft(x))^2)[1:floor((N-1)/2)]
-  
-  hprod <- function(x,y) sum(x * hN * y)/N
-  
-  # B_i = sum[m<= i](A_im) f_m the orthogonal base signals
-  # x_m what remains of the signal at time step m
-  
-  S <- rep(0,2*nfreq)
-  Prod <- rep(0,2*nfreq)
-  amp <- rep(NA,2*nfreq)
-  A <- matrix(0,2*nfreq,2*nfreq)
-  Q <- matrix(0,2*nfreq,2*nfreq)
-  f <-  list()
-  x <-  list()
-  B <- list()
-  # S[m] = hprod (f[m], B[m])
-  freqs = 2.*pi*seq(0,(N-1))/N
-  x[[1]] = xdata
-  
-   
-  for (m in seq(2*nfreq)){
-  hx <- hN*x[[m]]
-    # there is a little tweak here: 
-    # if we have reached 2*nfreq and not nu[m] is na, then
-    # it means that we have encountered a constant somewhere, and that last step should be skipped. 
-    if ( ! ( m == 2*nfreq && is.na(nu[m]))) {
-    # ok we are on business
-    if (is.na(nu[m])){
-    # is the "m" frequency already provided ? 
-    # either there were provided from the start and developped  in the first lines of this routine
-    # or either they were not provided, but we are in this case where it was set in "m-1" because
-    # we identified a non-null frequency
-    # in both configurations, we look at a frequency with the fourier transform
-    fbase <- freqs[which.max(power(hx))]
-    brackets <- c(fbase-pi/N, fbase+pi/N);
-    thx <- t(hx)
-    
-    # after profiling, the fastest seems the first option below
-    # this is the weak link
-    # coding the function in c might be an option
-    
-    fmax = cmna::goldsectmax(function(f) 
-                             {
-                              ft <- f*t
-                             (thx %*% cos(f*t))^2 + (thx %*% sin(f*t))^2}, 
-                             brackets[1], brackets[2], tol=1.e-10, m=9999)
-    
-   #  fmax = cmna::goldsectmax(function(f) {
-   #                           ft <- f*t
-   #                           (sum(hx * cos(ft)))^2 + (sum(hx %*% sin(ft)))^2}, 
-   #                           brackets[1], brackets[2], tol=1.e-10, m=9999)
-
-   #  fmax = cmna::goldsectmax(function(f) {
-   #                           ft <- f*t
-   #                           Mod(sum(hx * exp(1i*ft)))}, 
-   #                           brackets[1], brackets[2], tol=1.e-10, m=9999)
-
-
-
-    if (fmax > freqs[2]/2){
-      # if we really identified a frequency
-      # we are in this case where the frequency was not a priori provided
-      # we se set this phase to zero, and the next one to pi/2, and also set the frequencies accordinly
-    phase[m] <- 0.
-    nu[m]  <- fmax
-    phase[m+1] <- pi/2. 
-    nu[m+1] <- -fmax
-    } else {
-      # again we are still in thes case where a freqency was not a priori provided
-      # but the more particular case where the frequency is (with tolerance) considered as zero
-      # f[[m]] will effectively be constant. 
-      phase[m] <- 0. 
-      nu[m] <- 0.
-    }}
-    
-    f[[m]] <- cos(nu[m]*t + phase[m])
-
-    if (fast)  # we use analytical forms of the products
-    {
-       
-    for (i in seq(m))
-    {
-
-        num <- (nu[m] - nu[i])*T
-        nup <- (nu[m] + nu[i])*T
-        phim <- (phase[m] - phase[i])
-        phip <- (phase[m] + phase[i])
-
-        Qm <- ifelse(num == 0, 1 ,  Q(num))
-        Qp <- ifelse(nup == 0, 1 ,  Q(nup))
-        cosm <- cos(phim) * cos(num) * Qm
-        sinm <- sin(phim) * sin(num) * Qm
-        cosp <- cos(phip) * cos(nup) * Qp
-        sinp <- sin(phip) * sin(nup) * Qp
-
-        Q[m,i] = 0.5 * (cosm + cosp - sinm - sinp)
-    }
-    } else {
-    for (i in seq(m)) 
-    {Q[m,i] = hprod(f[[m]],f[[i]])  # so remember the convetion
-                                                 # these are symmetric matrices
-                                                 # and we only populate the j <= m 
-    }
-    }
-
-    #
-    # before normalisation
-    # B[[m]] = sum_j=1,m A[m,j]*f[j] et
-    # B[[m]] = ( f[m] - sum_(1,m-1) (f[m] %*% B[i]) * B[i]
-    # B[[m]] = ( f[m] - sum_(1,m-1) 
-    #                (f[m] * sum(j=1,i A_j,i  f[i]) ) *
-    #                               ( sum (A_j=1, i) A_j,i f[i] )
-    
-    #/ ||B[[m]]||
-    
-    # one can then verify that:
-    
-    # B[[m]] %*% B[j]] = 
-    # ( f[m] - sum_(1,m-1) (f[m] %*% B[i]) * B[i] ) %*%  B[j] =
-    # ( f[m] - (f[m] %*% B[j]) * B[j] ) %*%  B[j] =
-    # ( f[m] %*% B[j] - f[m] %*% B[j]) = 0 
-    
-    # before normalisation
-    A[m,] = 0
-    A[m, m] = 1.
-    if (m>1){
-      # fmbi = f[m] %*% B[i]
-      fmbi <- rep(0,(m-1))
-      for (j in seq(m-1)) for (i in seq(j)) fmbi[j] = fmbi[j] + A[j,i]*Q[m,i]
-      for (j in seq(m-1)) for (i in seq(j,(m-1))) A[m,j] = A[m,j] - fmbi[i]*A[i,j]
-    }
-    # so, at this stage, we have B[[m]] = sum [A[m,j]] * f[j]
-    
-    # B[[2]] = (A[2,1]*f[1] + A[2,2]*f[2])
-    # B[[3]] = (A[3,1]*f[1] + A[3,2]*f[2] + A[3,3] * f[3] ) 
-    
-
-    norm = 0
-    if (m > 1){
-    for (i in seq(m)) {
-      norm = norm + (A[m,i]*A[m,i]*Q[i,i])
-      if (i>1) for (j in seq(i-1)) { norm  = norm + 2*A[m,i]*A[m,j]*Q[i,j]
-    }}
-    } else {
-      norm <- A[m,m]*A[m,m]*Q[m,m]
-    }
- 
- 
-    
-    A[m,] = A[m,] / sqrt(norm)
-      
-     
-    Prod[m] = hprod(x[[1]],f[[m]])
-    # les prods sont corrects
-
-    S[m] = 0. 
-    for (j in seq(m))  S[m] = S[m] + Prod[j]*A[m,j]
-    # for (j in seq(m))  S[m] = S[m] + A[m,j] * hprod(x[[1]], f[[j]])
-
-    # les Sm sont les projections dans la nouvelle base
-    # not necessary, for verification only
-    # computes the B and verify they are orthonormal
-     B[[m]]=0
-     for (j in seq(m)) B[[m]] = B[[m]] + A[m,j] * f[[j]]
-
-
-    #print ('S[m] calcule de deux facons')
-    #print (m)
-    #print (S[m])
-    #print (hprod(x[[1]], B[[m]]))
-    #print ('end test')
-
-    # for (j in seq(m)) print (sprintf("B%i%i : %3g", i, j, hprod(B[[j]], B[[m]])))
-    # if you are curious, the amplitude of the signal will be
-    # amp[j] = contribution of the sum( S[m]*B[m]) to the f[[m]]
-    # amp[j] = sum(m in seq(j)) S(m) * A[m,j]
-    
-    # for (i in seq(m)) amp[j] = amp[j]+S[m]*A[m,j]
-    # and the residual signal
-    
-    # x[[m+1]] = x[[m]] - S[m] * B[[m]]
-
-
-    x[[m+1]] = x[[m]]
-    # for (j in seq(m))  x[[m+1]] = x[[m+1]] - S[m] * A[m,j]*f[[j]]
-    x[[m+1]] = x[[m+1]] - S[m] * B[[m]]
-#     print('residu final')
-#     print(m)
-#     print (x[[m+1]][seq(20)])
-  
-  }
-
-  }
-  
-  #print ('TESTESTEST')
-  #xtest <- 0 
-  #xtest2 <- 0 
-  #coeftests <- rep(0,m)
-  #phi1 <- 0.4 ; phi2 <- 0.1234
-  # coefreals <- c(cos(phi1), -sin(phi1), cos(phi2), -sin(phi2))
-  # for (j in seq(m)) xtest <- xtest + coefreals[j] * f[[j]]
-  # for (j in seq(m)) for (i in seq(j)) xtest2 <- xtest2 + S[j] * A[j,i] * f[[i]]
-  # for (j in seq(m)) for (i in seq(j)) coeftests[i] <- coeftests[i] + S[j] * A[j,i]
-  # print ((xtest - x[[1]])[seq(20)])
-  # print ('ENDTESTESTET1')
-  # print ((xtest2 - x[[1]])[seq(20)])
-  # print ('ENDTESTESTET2')
-  # print ("cofreals")
-  # print (coefreals)
-  # print ("coftests")
-  # print (coeftests)
-  
-  # donc: j'ai demontre que les deux xtests sont les memes
-  # coefreals est la bonne solution
-  # et je sais ussi que B = sum A_mj f_j
-  # donc je devroais simplement avour les Sj A_mj
-  
-
-  # at htis point I would like to make a littel check
-  # make A full
-  # is A Q t(A) = BtB orthonormal ? 
-  #print ('TEST')
-  #for (m in seq(2,2*nfreq)) for (j in seq(m-1)) A[j,m]=A[m,j]
-  #for (m in seq(2,2*nfreq)) for (j in seq(m-1)) Q[j,m]=Q[m,j]
-  #print (A*t(Q)*t(A))
-  print ('END TEST')
-
-  print (A) 
-  print (A %*% S)
-  # amplitudes
-  
-  mmax <- 2*nfreq
-  if (is.na(nu[mmax])) mmax <- mmax - 1
-
-  amp[1:mmax] = 0;
-  for (m in seq(mmax)) for (j in seq(m)) amp[j] <- amp[j] + S[m] * A[m,j]
-#   print ('ampp...')
-#   print (amp)
-
-  amp2m <- amp;
-
-  for (m in seq(mmax)){
-    # if the perivous "m" was already of the same frequency, we can merge amplitudes and set phase
-    if ((m > 1) && (nu[m-1] == -nu[m])){
-        
-        phase[m] <- Arg (-1i * amp[m] + amp[m-1])
-        amp[m] <- sqrt(amp[m-1]^2 + amp[m]^2)
-        print ("phase")
-        print (c(m, phase[m], amp[m], amp[m-1]))
-        amp[m-1] <- NA
-        phase[m-1] <- NA
-    }
-  }
-#   print ('amp2m')
-#   print (amp2m)
-#   print ('S')
-#   print (S)
-#   print ('PRod')
-#   print (Prod)
-  
-  # at this point we should have exactly nfreq non na values
-  # we print a message if this is not right, and in that case I suspect some debugging will be needed. 
-
-  valid_frequencies <- which(!is.na(amp))
-  nu <- -nu[valid_frequencies]
-  amp <- amp[valid_frequencies]
-  phase <- phase[valid_frequencies]
-  if (length(valid_frequencies) != nfreq) message (sprintf("something goes wrong : %i valid frequencies, and nfreq = %i", valid_frequencies, nfreq))
-#   print ("nu")
-#   print(nu)
-#   print ("amp")
-#   print (amp)
-
-  OUT = data.frame(nu=nu, amp=amp, phase=phase) 
-  return(OUT)
-}
-  
-     
-#' Modified Fourier transform  for real series (variant)
-#'
-#' R-coded version of the Modified Fourier Transform
-#' with frequency  correction, adapted to R. 
-#' much slower than mfft (for complex numbers) as the latter is
-#' mainly written in C, but is physically
-#' more interpretable if signal is real, because
-#' it is designed to have no imaginary part in the residual
-#' A C-version should be supplied one day. 
-#'
-#' @importFrom cmna goldsectmax
-#' @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.
-#' @param fast (default = TRUE) uses analytical formulations for the crossproducts involving sines and cosines
-#' @param correction:  0: no frequency correction (equivalent to Laskar); 1 : frequency correction using linear approximation ; 2: frequency correction using sythetic data
-#' @author Michel Crucifix
-#' @references
-#' \insertRef{sidlichovsky97aa}{gtseries}
-#' @examples
-#'
-#' set.seed(12413)
-#' t = seq(1024)
-#' x_orig = cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498) + rnorm(1024)*0.12
-#' OUT <- mfft_real(x_orig)
-#' print(OUT)
-#'
-#' @export mfft_real_quatro
-  # will withold the definitive frequencies
-mfft_real_quatro <- function(xdata, nfreq=5, correction=1, fast=TRUE){
-    N <- length(xdata)
-    N2 <- N/2.
-    xdata = stats::as.ts(xdata)
-    dt = deltat(xdata)
-    startx = stats::start(xdata)[1]
-    N <- length(xdata)
-    OUT <- analyse_quatro(xdata, nfreq, fast)
-
-
-    # correction  (methode 2)
-#   }
-  if (correction == 2){
-   xdata_synthetic <- rep(0,N)
-   t <- seq(N)-1
-   for (i in seq(nfreq)) xdata_synthetic = xdata_synthetic + OUT$amp[i]*cos(OUT$nu[i]*t + OUT$phase[i])
-   OUT2 <- analyse_quatro(xdata_synthetic, nfreq, fast)
-   OUT$nu = OUT$nu + (OUT$nu - OUT2$nu)
-   OUT$amp = OUT$amp + (OUT$amp - OUT2$amp)
-   OUT$phase = OUT$phase + (OUT$phase - OUT2$phase)
-   } else if (correction == 1){
-
-    for (j in seq(nfreq)){
-       epsilon = OUT$amp[j] * Qprime(-2 * OUT$nu[j] * N2)*cos(2 * OUT$nu[j] * N2 + 2 * OUT$phase[j])
-#        print ('epsilon')
-#        print (c(epsilon,j, OUT$amp[j], OUT$nu[j], OUT$phase[j]))
-       if ((j+1) <= nfreq) { for (s in seq(j+1, nfreq)) {
-        epsilon = epsilon + OUT$amp[s]  * 
-         ( 
-          Qprime( (OUT$nu[s] - OUT$nu[j])*N2)*cos((OUT$nu[j] - OUT$nu[s])*N2 + OUT$phase[j] - OUT$phase[s] ) -
-          Qprime(( OUT$nu[s] + OUT$nu[j])*N2)*cos((OUT$nu[j] + OUT$nu[s])*N2 + OUT$phase[j] + OUT$phase[s] ) )
-     }}
-#     print(sprintf("j= %i, OUT$nu[j]= %.4f, OUT$amp[j] = %.4f", j, OUT$nu[j], OUT$amp[j]))   
-    epsilon = epsilon / Qsecond0 / N2 / OUT$amp[j]
-
-    OUT$nu[j] = OUT$nu[j] - epsilon
-    }
-    OUT <- analyse_quatro(xdata, nfreq, fast, nu = OUT$nu)
-  }
-
-  # account for tseries scaling (i.e. uses the vaue of deltat and start encoded
-  # in the time series object 
-    OUT$nu <- OUT$nu / dt
-    OUT$phase <- OUT$phase - startx*OUT$nu
-
-  # rearrange terms to avoid negative amplitudes and have phases in the right quandrant
-  # these are cosines so even functions 
-
-    to_be_corrected <- which (OUT$amp < 0)
-    if (length(to_be_corrected)){
-      OUT$amp[to_be_corrected] <- - OUT$amp[to_be_corrected]
-      OUT$phase[to_be_corrected] <- OUT$phase[to_be_corrected] + pi 
-    }
-
-    to_be_corrected <- which (OUT$nu < 0)
-    if (length(to_be_corrected)){
-      OUT$nu[to_be_corrected] <- - OUT$nu[to_be_corrected]
-      OUT$phase[to_be_corrected] <- - OUT$phase[to_be_corrected] 
-    }
-
-   # finally, order termes by decreasing amplitude 
-    o <- order(OUT$amp, decreasing = TRUE)
-    OUT$amp <- OUT$amp[o]
-    OUT$nu <- OUT$nu[o]
-    OUT$phase <- OUT$phase[o]
-
-    OUT$phase <- (OUT$phase + (2*pi)) %% (2*pi)
-
-
-    # rename for class compatibility
-    names(OUT) <- c("Freq","Amp","Phases")
-    attr(OUT, "class") <- "mfft_deco"
-    attr(OUT, "data")  <- xdata
-    attr(OUT, "nfreq")  <- nfreq
-    return(OUT)
-}
-

TODO

-- make nicer graphs for mfft  <- in prgoress
-- finalise computations for mfft_real_ter (use notebook) <- seems to be ok forfrequenices
-- there still seems to be a bias on amplitudes that is not well corrected with the analytical methods, when there are several frequencies. As if the first components where eating a bit too much power. Not sure why, and not sure how to correct this. 
- But the bias is much less important with the complex version. So there must be a good mathematical reason that is not
- obvious to me (in mfft_complex not quite the same power goes in omega and -omega, so maybe the explanation is there ?)
-- use that one only; deprecate mfft_real
+- make nicer graphs for mfft  <- in progress
 - work on documentation
+- document data and add astronomical data (documentation has to be revised for harmonic and pseudo_log)
+- think more about the significance spectrum
+- revis the haar bit
+- revise the periodogram
 - and make a paper ! 

deprecated/mfft_real.R

+calc_amp <- function(x, nfreq, Fs){
+  N <- length(x)
+  t <- seq(N)
+  xx = matrix(rep(1,N),N,1)
+  freqs = 2.*pi*seq(0,(N-1))/N
+  power <- function(x) (Mod(fft(x))^2)[1:floor((N-1)/2)]
+  residu <- x  - mean(x)
+  for (j in seq(nfreq)) {
+ 
+# temps de calcul peut etre gagne en C
+
+
+    if (Fs[j] == (-1) ){
+       hresidu = gsignal::hann(N)*residu
+       fbase <- freqs[which.max(power(hresidu))]
+       brackets <- c(fbase-pi/N, fbase+pi/N);
+       thresidu <- t(hresidu)
+       fmax = cmna::goldsectmax(function(f) 
+                          (thresidu %*% cos(f*t))^2 + (thresidu %*% sin(f*t))^2, 
+                          brackets[1], brackets[2], tol=1.e-10, m=9999)
+       Fs[j] = fmax
+    } 
+ 
+     xx <- cbind(xx, cos(Fs[j]*t), sin(Fs[j]*t))
+     
+# cette partie doit pouvoir etre acceleree en tirant partit du fait
+# que l'expression analytitque de xx%*xx est connue
+# voir non-working code fmft_real.C
+# Il faut utiliser l'integrale connue int_0^2pi cos(x)hann(x) 
+# on utilise: cos(a)cos(b) = 0.5*(cos(a+b) - cos(a-b)) etc pour toutes les comb. cossin
+# et int_0^2T cos(wx) = 1/W(sin[2T]-sin[0]) , ou si on applique le produit scalaire avec fenetre Hann
+# et int_0^2T cos(wx)Hann[x] = (pi^2)(pi^2-w^2T^2)(sin(wT)/wT) * (coswT)
+# et int_0^2T sin(wx)Hann[x] = (pi^2)(pi^2-w^2T^2)(sin(wT)/wT) * (sin(wT))
+ # sin(wT)/(wT) = 1
+ # et on se rappelle que la decomposition gramschmidt
+ # c'est simplement, avec x_i le ie vecteur
+ # xn_1 = x_1 / ||x_1||  
+ # xn_(i+1) = x_(i+1) - sum_(j=1,i) x_(i+1)%*%x(i)/||x_(i+1)||
+# en C, cela doit reduire le temps de calcul considerablement
+
+
+     B <- pracma::gramSchmidt(xx)$Q
+     Coefs = t(xx) %*% B
+     aa <- as.numeric(residu) %*% B
+
+     signal <- B %*% t(aa)
+     residu <- residu - signal
+  }
+  
+  aa <- as.numeric(x) %*% B
+  sol = solve(t(Coefs), t(aa))
+  sol1 <- matrix(sol[-1],2, nfreq)
+  
+  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(Freq=Freqs, Amp=amps, Phases=Phases) 
+  attr(OUT,"class") = "mfft_deco"
+  return(OUT)
+}
+
+
+
+#' Modified Fourier transform  for real series
+#'
+#' R-coded version of the Modified Fourier Transform
+#' with frequency  correction, adapted to R. 
+#' much slower than mfft (for complex numbers) as the latter is
+#' mainly written in C, but is physically
+#' more interpretable if signal is real, because
+#' it is designed to have no imaginary part in the residual
+#' A C-version should be supplied one day. 
+#'
+#' @importFrom cmna goldsectmax
+#' @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.
+#' @author Michel Crucifix
+#' @references
+#' \insertRef{sidlichovsky97aa}{gtseries}
+#' @examples
+#'
+#' set.seed(12413)
+#' t = seq(1024)
+#' x_orig = cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498) + rnorm(1024)*0.12
+#' OUT <- mfft_real(x_orig)
+#' print(OUT)
+#'
+#' @export mfft_real
+mfft_real <- function(xdata, nfreq=5, correction=TRUE){
+  xdata = stats::as.ts(xdata)
+  dt = deltat(xdata)
+  startx = stats::start(xdata)[1]
+  N = length(xdata)
+  times <- (seq(N)-1)*dt+startx
+  residu = xdata
+
+  # will withold the definitive frequencies
+  Fs <- rep(-1, nfreq)
+    k=1
+    OUT <- calc_amp(xdata, nfreq, Fs)
+
+    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) {
+##         (1 / y) * (pi^2/ (pi^2-y^2))*(cos(y) + sin(y) / y * (3. * y^2 - pi^2)/(pi^2 - y^2))}
+##       Q0 <- {2. / (pi^2) -1./3.}
+##       
+##       epsilon = rep(0,nfreq)
+##       
+## #       tau = (N-1) / 2.
+##      
+##       for (j in seq(k,nfreq-1)){
+##          jj = j+1
+##          ysj =(Freqs[(jj+1):(nfreq+1)]-Freqs[jj])*tau
+##          epsilon[j] = sum(amps[(jj + 1):(nfreq+1)]*Qp(ysj)/
+##                           (amps[jj] * Q0 * tau) * cos(ysj + Phases[(jj+1):(nfreq+1)]-Phases[jj]))
+##       }
+##       Epsilon=c(0,epsilon)
+##       Freqs <- Freqs - Epsilon
+##       Fs <- Freqs[2:(nfreq+1)]
+##     }
+## 
+
+   
+    # correction  (methode 2)
+
+#   }
+  if (correction){
+   xdata_synthetic <- rep(0,N)
+     Fs <- rep(-1, nfreq)
+    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$Freq - OUT2$Freq)
+    # adjust to units
+    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
+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))
+}
+
+
+

man/mfft.Rd

 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/mfft.R
+% Please edit documentation in R/mfft_complex.R
 \name{mfft}
 \alias{mfft}
 \title{Modified Fourier transform}
 \usage{
-mfft(xdata, minfreq = NULL, maxfreq = NULL, flag = 1, nfreq = 30)
+mfft(xdata, minfreq = NULL, maxfreq = NULL, correction = 1, nfreq = 30)
 }
 \arguments{
 \item{xdata}{The data provided either as a time series (advised), or as a vector. 
@@ -12,30 +12,40 @@ may be complex}
 
 \item{minfreq, maxfreq}{If provided, bracket the frequencies to be probed. Note this are
 more precisely angular velocities (2\pi / period), expressed in time-inverse units
-with the time resolution encoded in `xdata` if the latter is a time series.}
-
-\item{flag}{:
-     Modified Fourier Transform                  if   flag = 1;
-     Frequency Modified Fourier Transform        if   flag = 2;
-     FMFT with additional non-linear correction  if   flag = 3
-(while the first algorithm is app. 3 times faster than the third one, 
-the third algorithm should be in general much more precise).  
+with the time resolution encoded in `xdata` if the latter is a time series. 
 The computed frequencies are in the range given by minfreq and maxfreq.}
 
+\item{correction}{:
+The parameter is passed on to `mfft_real` or `mfft_complex` depending on
+`xdata` is real or complex, but the meaning of the parameter is not
+quite the same for both. This needs to be fixed. 
+one approach would be to define scalars like 'fft', 'mfft', 'mfft_linear', 
+'mfft_second_order_correction'.}
+
 \item{nfreq}{is the number of frequencies returned, must be smaller that the length of  xdata.}
 }
 \value{
-STILL NEED TO BE DESCRIBED
+a `mfft_deco` object, based on a data.frame with columns "Freq", "Ampl" and "Phases". 
+        note that because of a language glitch (to be fixed), "Freq" actually means "Rate"
 }
 \description{
 Implementation of the the Frequency Modified Fourier Transform 
 (Sidlichovsky and Nesvorny 1997, Cel. Mech. 65, 137). 
-Given a quasi--periodic complex signal X + iY, the algorithm 
+Given a quasi-periodic complex signal X + iY, the algorithm 
 estimates the frequencies (f_j), amplitudes (A_j) and phases 
-(psi_j) in its decomposition.
+(psi_j) in its decomposition. The algorithm was generalised
+by M. Crucifix when the input signal is real. 
+`mfft` is a wrapper thar redirects the data
+to `mfft_real` or `mfft_complex` depending on the nature of
+the input data. 
+TODO: this needs more work  because the two routines
+do not yet take exactly the same inputs. In particular, the meaning
+the `correction` parameter is not the same 
+TODO: add Laskar reference and say that the correction=0 version
+      corresponds to Laskar's routine.
 }
 \author{
 Michel Crucifix for the R code, and David Nesvorny for most of the supporting C code doing the
-actual computations
+        Frequency Modified Fourier transform  for Complex Numbers
 \insertRef{sidlichovsky97aa}{gtseries}
 }

man/mfft_deco.Rd

 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/mfft_real.R
+% Please edit documentation in R/mfft_support.R
 \name{reconstruct_mfft}
 \alias{reconstruct_mfft}
 \alias{mfft_anova}
@@ -25,4 +25,5 @@ mfft_anova(M)
 MFFT reconstruction
 
 MFFT ANOVA
+not ready. do not use.
 }

man/mfft_real.Rd

@@ -2,14 +2,33 @@
 % Please edit documentation in R/mfft_real.R
 \name{mfft_real}
 \alias{mfft_real}
-\title{Modified Fourier transform  for real series}
+\title{Frequency Modified Fourier transform  for real series}
 \usage{
-mfft_real(xdata, nfreq = 5, correction = TRUE)
+mfft_real(
+  xdata,
+  nfreq = 5,
+  minfreq = NULL,
+  maxfreq = NULL,
+  correction = 1,
+  fast = TRUE
+)
 }
 \arguments{
 \item{xdata}{The data provided either as a time series (advised), or as a vector.}
 
 \item{nfreq}{is the number of frequencies returned, must be smaller that the length of  xdata.}
+
+\item{minfreq, maxfreq}{If provided, bracket the frequencies to be probed. Note this are
+more precisely angular velocities (2\pi / period), expressed in time-inverse units
+with the time resolution encoded in `xdata` if the latter is a time series.}
+
+\item{fast}{(default = TRUE) uses analytical formulations for the crossproducts involving sines and cosines. 
+note: this is not really faster because the bottleneck is actually the goden section search. But more elegant.}
+
+\item{correction:}{0: no frequency correction (equivalent to Laskar); 1 : frequency correction using linear approximation ; 2: frequency correction using sythetic data (both are documented in the Sidlichovsky and Nesvorny paper.}
+}
+\value{
+a `mfft_deco` object, based on a data.frame with columns "Freq", "Ampl" and "Phases".
 }
 \description{
 R-coded version of the Modified Fourier Transform
@@ -22,11 +41,9 @@ A C-version should be supplied one day.
 }
 \examples{
 
-set.seed(12413)
-t = seq(1024)
-x_orig = cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498) + rnorm(1024)*0.12
-OUT <- mfft_real(x_orig)
-print(OUT)
+data(harmonic_sample)
+spectrum <- mfft_real(harmonic_sample$data)
+print(spectrum)
 
 }
 \references{

man/mfft_real_C.Rd

-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/mfft_real_C.R
-\name{mfft_real_C}
-\alias{mfft_real_C}
-\title{Modified Fourier transform (Real C)}
-\usage{
-mfft_real_C(xdata, minfreq = NULL, maxfreq = NULL, flag = 1, nfreq = 30)
-}
-\arguments{
-\item{xdata}{The data provided either as a time series (advised), or as a vector. 
-may be complex}
-
-\item{minfreq, maxfreq}{If provided, bracket the frequencies to be probed. Note this are
-more precisely angular velocities (2\pi / period), expressed in time-inverse units
-with the time resolution encoded in `xdata` if the latter is a time series.}
-
-\item{flag}{:
-     Modified Fourier Transform                  if   flag = 1;
-     Frequency Modified Fourier Transform        if   flag = 2;
-     FMFT with additional non-linear correction  if   flag = 3
-(while the first algorithm is app. 3 times faster than the third one, 
-the third algorithm should be in general much more precise).  
-The computed frequencies are in the range given by minfreq and maxfreq.}
-
-\item{nfreq}{is the number of frequencies returned, must be smaller that the length of  xdata.}
-}
-\value{
-STILL NEED TO BE DESCRIBED
-}
-\description{
-Implementation of the the Frequency Modified Fourier Transform 
-(Sidlichovsky and Nesvorny 1997, Cel. Mech. 65, 137). 
-Given a quasi--periodic complex signal X + iY, the algorithm 
-estimates the frequencies (f_j), amplitudes (A_j) and phases 
-(psi_j) in its decomposition.
-}
-\author{
-Michel Crucifix for the R code, and David Nesvorny for most of the supporting C code doing the
-actual computations
-\insertRef{sidlichovsky97aa}{gtseries}
-}

man/mfft_real_quatro.Rd

-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/mfft_real_quatro.R
-\name{mfft_real_quatro}
-\alias{mfft_real_quatro}
-\title{Modified Fourier transform  for real series (variant)}
-\usage{
-mfft_real_quatro(xdata, nfreq = 5, correction = 1, fast = TRUE)
-}
-\arguments{
-\item{xdata}{The data provided either as a time series (advised), or as a vector.}
-
-\item{nfreq}{is the number of frequencies returned, must be smaller that the length of  xdata.}
-
-\item{fast}{(default = TRUE) uses analytical formulations for the crossproducts involving sines and cosines}
-
-\item{correction:}{0: no frequency correction (equivalent to Laskar); 1 : frequency correction using linear approximation ; 2: frequency correction using sythetic data}
-}
-\description{
-R-coded version of the Modified Fourier Transform
-with frequency  correction, adapted to R. 
-much slower than mfft (for complex numbers) as the latter is
-mainly written in C, but is physically
-more interpretable if signal is real, because
-it is designed to have no imaginary part in the residual
-A C-version should be supplied one day.
-}
-\examples{
-
-set.seed(12413)
-t = seq(1024)
-x_orig = cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498) + rnorm(1024)*0.12
-OUT <- mfft_real(x_orig)
-print(OUT)
-
-}
-\references{
-\insertRef{sidlichovsky97aa}{gtseries}
-}
-\author{
-Michel Crucifix
-}

man/mfft_real_ter.Rd

-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/mfft_real_ter.R
-\name{mfft_real_ter}
-\alias{mfft_real_ter}
-\title{Modified Fourier transform  for real series (variant)}
-\usage{
-mfft_real_ter(xdata, nfreq = 5, correction = 1, fast = TRUE)
-}
-\arguments{
-\item{xdata}{The data provided either as a time series (advised), or as a vector.}
-
-\item{nfreq}{is the number of frequencies returned, must be smaller that the length of  xdata.}
-
-\item{fast}{(default = TRUE) uses analytical formulations for the crossproducts involving sines and cosines}
-
-\item{correction:}{0: no frequency correction (equivalent to Laskar); 1 : frequency correction using linear approximation ; 2: frequency correction using sythetic data}
-}
-\description{
-R-coded version of the Modified Fourier Transform
-with frequency  correction, adapted to R. 
-much slower than mfft (for complex numbers) as the latter is
-mainly written in C, but is physically
-more interpretable if signal is real, because
-it is designed to have no imaginary part in the residual
-A C-version should be supplied one day.
-}
-\examples{
-
-set.seed(12413)
-t = seq(1024)
-x_orig = cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498) + rnorm(1024)*0.12
-OUT <- mfft_real(x_orig)
-print(OUT)
-
-}
-\references{
-\insertRef{sidlichovsky97aa}{gtseries}
-}
-\author{
-Michel Crucifix
-}