Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
#' 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")
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))
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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
#
#
# 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)