Skip to content
Extraits de code Groupes Projets
Valider 79cb39de rédigé par Michel Crucifix's avatar Michel Crucifix
Parcourir les fichiers

add mmfft

parent 63767df7
Aucune branche associée trouvée
Aucune étiquette associée trouvée
Aucune requête de fusion associée trouvée
......@@ -14,9 +14,11 @@ export(attributeTones)
export(cwt_morlet)
export(hilbert_extension)
export(mem)
export(mfft)
export(mfft_anova)
export(mfft_complex)
export(mfft_real)
export(mmfft)
export(periodogram)
export(powerspectrum.wavelet)
export(reconstruct_mfft)
......
......@@ -9,9 +9,6 @@
#' `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.
......@@ -20,20 +17,19 @@
#' 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 correction: 0: no frequency correction (equivalent to Laskar); 1 : frequency correction using linear approximation ; 2: frequency correction using sythetic data;
#' 3: second order-correction using synthetic data (all documented in the Sidlichovsky and Nesvorny reference). Note: 1 is not implemented for complex time series; 4 is
#' not documented for real time series
#' @param nfreq is the number of frequencies returned, must be smaller that the length of xdata.
#' @param force_complex : use the complex number implementation even if the time series is real.
#' @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
#' @export mfft
mfft <- function(xdata, nfreq=15, minfreq=NULL, maxfreq=NULL, correction=1, force_complex = FALSE) {
if (is.complex(xdata) || force_complex) return(mfft_complex(xdata, nfreq, minfreq, maxfreq, correction)) else
return(mfft_real(xdata, nfreq, minfreq, maxfreq, correction))
}
......@@ -76,10 +72,8 @@ mfft_complex <- function(xdata, nfreq=30, minfreq=NULL, maxfreq=NULL, correctio
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")
flag <- c(1,'NA',2,3)[correction+1]
if (is.na(flag)) stop('this correction scheme is not implemented for complex time series')
if (is.null(minfreq)){
my_minfreq <- -pi
......@@ -108,21 +102,15 @@ mfft_complex <- function(xdata, nfreq=30, minfreq=NULL, maxfreq=NULL, correctio
# 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"
class(OUT) <- c("mfft_deco", "data.frame")
attr(OUT,"nfreq") <- nfreq
attr(OUT,"data") <- xdata
return(OUT)
}
#
#
# return(OUTVEC)
# dyn.load('fmft.so')
......
......@@ -337,7 +337,8 @@ mfft_analyse <- function(xdata, nfreq, fast = TRUE, nu = NULL, minfreq=NULL, max
#' @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: 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.
#' @param correction: 0: no frequency correction (equivalent to Laskar); 1 : frequency correction using linear approximation ; 2: frequency correction using sythetic data;
#' 3: second order-correction using synthetic data (all documented in the Sidlichovsky and Nesvorny reference)
#' @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.
#' note: this is not really faster because the bottleneck is actually the goden section search. But more elegant.
......@@ -354,6 +355,8 @@ mfft_analyse <- function(xdata, nfreq, fast = TRUE, nu = NULL, minfreq=NULL, max
#' @export mfft_real
mfft_real <- function(xdata, nfreq=5, minfreq=NULL, maxfreq=NULL, correction = 1 , fast=TRUE){
if (correction == 4) "this correction scheme is currently not implemented for real time series"
N <- length(xdata)
N2 <- N/2.
xdata = stats::as.ts(xdata)
......@@ -427,7 +430,7 @@ mfft_real <- function(xdata, nfreq=5, minfreq=NULL, maxfreq=NULL, correction =
# rename for class compatibility
names(OUT) <- c("Freq","Amp","Phases")
attr(OUT, "class") <- "mfft_deco"
class(OUT) <- c("mfft_deco", "data.frame")
attr(OUT, "data") <- xdata
attr(OUT, "nfreq") <- nfreq
return(OUT)
......
......@@ -52,9 +52,9 @@ plot.mfft_deco <- function (M,periods=FALSE,labels=NULL,...){
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)
plabels <- as.character(1/frequencies)
if (0 %in% frequencies) plabels[which(frequencies == 0)] = "∞"
axis(1, line=3, at=2*pi*frequencies, labels=plabels)
mtext("Rate", 1, 2)
mtext("Period", 1, 4)
} else {
......@@ -62,8 +62,8 @@ plot.mfft_deco <- function (M,periods=FALSE,labels=NULL,...){
}
points(abs(M$Freq), abs(M$Amp),'p',...)
if (!is.null(labels)) {
yshift <- 0.05*range(M$Amp)
text(M$Freq, M$Amp + yshift, labels, srt=90, pos=4)
yshift <- 0.05*diff(range(M$Amp))
text(M$Freq, M$Amp, labels, srt=90, , adj=-0.4)
}
}
......
#' Moving Modified Fourier Transform
#'
#' @param xdata The data provided either as a time series (advised), or as a vector.
#' may be complex
#' @param seglength Length of moving segment. Both xdata length and seglenthm are adviced to be a power of 2 (faster)
#' @param ... passed to mmft
#' @return a `mmfft` object
#' @note in the current implementation, the right-hand-side of xdata that is left after an integer division
#' by seglenth is discarded
#' @author Michel Crucifix
#' @export mmfft
mmfft <- function(xdata, seglength = length(xdata) %/% 16, ...){
N <- length(xdata)
nsections <- (N %/% seglength)
N0 <- nsections * seglength
# remainder <- seglength - N0
xdata0 <- matrix( xdata[1:N0], byrow=FALSE, seglength, nsections)
myMfft <- function(x) mfft(x, ...)
OUT <- apply(xdata0, 2, myMfft)
class(OUT) <- "mmfft"
attr(OUT, "nsections") <- nsections
attr(OUT, "seglength") <- seglength
attr(OUT, "start") <- ifelse(is.ts(xdata), stats::start(xdata), 1)
attr(OUT, "deltat") <- ifelse(is.ts(xdata), stats::deltat(xdata), 1)
return(OUT)
}
#' @rdname mmfft
plot.mmfft <- function(x){
freqrange <-
c(min(sapply(x, function(xs) min(xs$Freq))),
max(sapply(x, function(xs) max(xs$Freq))))
amprange <-
c(min(sapply(x, function(xs) min(xs$Amp))),
max(sapply(x, function(xs) max(xs$Amp))))
print(amprange)
amp2lwd <- function(amp){ 3*amp/amprange[2] }
length <- attr(x, "nsections") *
attr(x, "seglength")
nsec <- length(x)
dt <- attr(x,"deltat")
tsec <- attr(x,"seglength")*dt
tstart <- attr(x,"start")
tend <- tstart + nsec * tsec
print(freqrange)
print(c(tstart, tend))
plot(c(tstart, tend), freqrange, type='n', xlab='Time', ylab='Rate')
for (iseq in seq(nsec)){
trange <- tstart + c(iseq-1,iseq)*tsec
obj <- x[[iseq]]
nfreq <- length(obj$Freq)
print(obj$Amp)
lwds <- sapply(obj$Amp, amp2lwd)
print(lwds)
for (j in seq(nfreq)){
lines(trange, rep(obj$Freq[j],2), lwd=lwds[j])
}
}
}
......@@ -4,7 +4,14 @@
\alias{mfft}
\title{Modified Fourier transform}
\usage{
mfft(xdata, nfreq = 30, minfreq = NULL, maxfreq = NULL, correction = 1)
mfft(
xdata,
nfreq = 15,
minfreq = NULL,
maxfreq = NULL,
correction = 1,
force_complex = FALSE
)
}
\arguments{
\item{xdata}{The data provided either as a time series (advised), or as a vector.
......@@ -17,12 +24,11 @@ more precisely angular velocities (2\pi / period), expressed in time-inverse uni
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{force_complex}{: use the complex number implementation even if the time series is real.}
\item{correction:}{0: no frequency correction (equivalent to Laskar); 1 : frequency correction using linear approximation ; 2: frequency correction using sythetic data;
3: second order-correction using synthetic data (all documented in the Sidlichovsky and Nesvorny reference). Note: 1 is not implemented for complex time series; 4 is
not documented for real time series}
}
\value{
a `mfft_deco` object, based on a data.frame with columns "Freq", "Ampl" and "Phases".
......@@ -38,9 +44,6 @@ 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.
}
......
......@@ -15,7 +15,7 @@ mfft_anova(M)
\method{as.data.frame}{mfft_deco}(x)
\method{plot}{mfft_deco}(M, periods = FALSE, ...)
\method{plot}{mfft_deco}(M, periods = FALSE, labels = NULL, ...)
\method{lines}{mfft_deco}(M, ...)
......@@ -25,6 +25,12 @@ mfft_anova(M)
\item{M}{: mfft_deco object}
\item{sum}{: TRUE if user wants to sum components in the reconstruction}
\item{periods}{if TRUE will add a lower axis with period labels}
\item{labels}{to be set above the frequency peaks. Can be the output of `attributeTone`}
\item{a}{`mfft_deco` object, typically the output of a `mfft` call.}
}
\value{
list of reconstructed components if sum=FALSE, full
......
......@@ -25,7 +25,8 @@ 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.}
\item{correction:}{0: no frequency correction (equivalent to Laskar); 1 : frequency correction using linear approximation ; 2: frequency correction using sythetic data;
3: second order-correction using synthetic data (all documented in the Sidlichovsky and Nesvorny reference)}
}
\value{
a `mfft_deco` object, based on a data.frame with columns "Freq", "Ampl" and "Phases".
......
......@@ -21,7 +21,7 @@ explicit label names, up to order 3 (this could be made more flexible is the fut
}
\examples{
omegas <- c( 0.123, 0.14312, 0.33251, 0.554313)
print(toneCombinations(omegas)
print(toneCombinations(omegas))
}
\author{
Michel Crucifix
......
Ce diff est replié.
0% Chargement en cours ou .
You are about to add 0 people to the discussion. Proceed with caution.
Terminez d'abord l'édition de ce message.
Veuillez vous inscrire ou vous pour commenter