-
Michel Crucifix a rédigéMichel Crucifix a rédigé
mfft_complex.R 8,38 Kio
#' 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))
}
#' Frequency Modified Fourier transform for Complex Numbers
#'
#' 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.
#' @useDynLib gtseries
#' @importFrom stats start
#' @importFrom stats as.ts
#' @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.
#' @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 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_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);
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 (is.null(minfreq)){
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))
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),
as.integer(ndata), Ixdata, Iydata, signal1,signal2,signal3, DUP=TRUE)
OUTVEC <- t(matrix(OUT[[7+flag]], 3, nfreq))
Freq <- OUTVEC[seq(nfreq) ]/dt
Ampl <- OUTVEC[nfreq+seq(nfreq) ]
Phase <- OUTVEC[2*nfreq+seq(nfreq) ] - startx*Freq
# 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)
}
#
#
# return(OUTVEC)
# dyn.load('fmft.so')
# ndata=8192
#for (exp in c('d')){
# namein <- file.path('La2010',sprintf("La2010%s_alkhqp3L.dat",exp));
# nameout <- sprintf("La2010%s_out.RData",exp);
# A <- read.table(namein)
# deltat=1000
# colnames(A) <- c('times','a','l','k','h','q','p')
# # ndata = as.integer64(ndata)
#
# times <- seq(ndata)
## nfreq <- as.integer(30)
# flag = 2
#
# minfreq <- as.double ( -1e6 / 180. / 3600. * pi /5. ) ;
# maxfreq <- as.double ( 1.e6 / 180. / 3600. * pi /5. ) ;
#
# signal1 = as.double (rep(0, nfreq*3))
# signal2 = as.double (rep(0, nfreq*3))
# signal3 = as.double (rep(0, nfreq*3))
#
# ntrunks <- floor(length(A$times)/ndata)
#
# # EMFFT <- sapply (seq(0,floor(length(A$times)/ndata)-1), function(i)
#
#
# timemin <- sapply (seq(0,ntrunks-1), function(i) { A$times [ndata* i + 1] });
# timemax <- sapply (seq(0,ntrunks-1), function(i) { A$times [ndata* i + ndata ] });
#
# EMFFT <- sapply (seq(0,ntrunks-1), function(i)
# {
# xdata <- A$h[ndata*i + (1:ndata)]
# ydata <- A$k[ndata*i + (1:ndata)]
#
# OUT = .C("fmft", nfreq, minfreq, maxfreq, as.integer(flag), as.integer(ndata), xdata, ydata, signal1,signal2,signal3, DUP=FALSE)
#
# # OUTVEC <- matrix(OUT[[8]], 10,3)
# OUTVEC <- t(matrix(OUT[[8]], 3, nfreq))
#
# return(OUTVEC)
#
# })
#
#
# # repack
#
# # repack
#
#
# Freq <- EMFFT[seq(nfreq), ]*360*60*60/2/pi/deltat
# Ampl <- EMFFT[nfreq+seq(nfreq), ]
# Phas <- EMFFT[2*nfreq+seq(nfreq), ]
#
#
# snake <- function(x1, x2, tol=2.e-1, mynfreq=nfreq){
# xnew <- c();
# xamp <- c();
# xpha <- c();
# xf1 <- x1[seq(mynfreq)]
# xf2 <- x2[seq(mynfreq)]
# xa1 <- x1[mynfreq+seq(mynfreq)]
# xa2 <- x2[mynfreq+seq(mynfreq)]
# xp1 <- x1[2*mynfreq+seq(mynfreq)]
# xp2 <- x2[2*mynfreq+seq(mynfreq)]
#
# while(length(xf1) > 0) {
# i <- which.min(abs(xf2 - xf1[1])+0*abs(xa2-xa1[1]));
# xnew <- c(xnew, xf2[i])
# xamp <- c(xamp, xa2[i])
# xpha <- c(xpha, xp2[i])
# xf2 <- xf2[-i]
# xf1 <- xf1[-1]
# xa2 <- xa2[-i]
# xa1 <- xa1[-1]
# xp2 <- xp2[-i]
# xp1 <- xp1[-1]
# }
# return(c(xnew, xamp, xpha))
# }
#
# mylist <- apply(EMFFT,2,as.numeric, simplify=FALSE)#
# terms <- simplify2array(Reduce(snake, mylist, mylist[[1]], accumulate=TRUE))
#
# myfreq <- terms[seq(nfreq),-1] * 360*60*60/2/pi/deltat
# myamp <- terms[nfreq+seq(nfreq),-1]
# myphase <- terms[2*nfreq+seq(nfreq),-1]
#
#
# rownames (myfreq) <- timemin
# rownames (myamp) <- timemin
# rownames (myphase)<- timemin
#
#
# OUT = list(Freq=myfreq, Amp=myamp, Phase=myphase)
# save(OUT, file=nameout)
#}: