diff --git a/NAMESPACE b/NAMESPACE index d9c5d5310b8124799a9eb75246be19ac1cae602e..64c9449fab41f390b41ccbd8e09e9ec5a539bde9 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -16,7 +16,6 @@ export(mem) export(mfft_anova) export(mfft_complex) export(mfft_real) -export(mfft_real_C) export(periodogram) export(powerspectrum.wavelet) export(reconstruct_mfft) diff --git a/R/mfft_complex.R b/R/mfft_complex.R index 9d9d5cc8e5b57bc3f96afb67204027844f51eb03..366d860c1a137fa8b560f5c47b2d408563ee1ef7 100644 --- a/R/mfft_complex.R +++ b/R/mfft_complex.R @@ -113,7 +113,7 @@ mfft_complex <- function(xdata, nfreq=30, minfreq=NULL, maxfreq=NULL, correctio attr(OUT,"class") <- "mfft_deco" attr(OUT,"nfreq") <- nfreq - attr(OUT,"xdata") <- xdata + attr(OUT,"data") <- xdata return(OUT) } diff --git a/R/mfft_support.R b/R/mfft_support.R index 22df4aee54bf0719a5f7f968630d75b263718ec4..b830019d1309d8e00dd43f086a42ea13ff2a4d2e 100644 --- a/R/mfft_support.R +++ b/R/mfft_support.R @@ -1,12 +1,19 @@ #' MFFT reconstruction +#' @importFrom stats deltat start #' @rdname mfft_deco +#' @param M : mfft_deco object +#' @param sum : TRUE if user wants to sum components in the reconstruction #' @export reconstruct_mfft -reconstruct_mfft <- function(M){ +#' @return list of reconstructed components if sum=FALSE, full +#' reconstructed time series otherwise +reconstruct_mfft <- function(M, sum=TRUE){ if (!(attr(M,"class")=="mfft_deco")) stop ("object is not a MFFT decomposition") - xdata <- attr(M,"xdata") + xdata <- attr(M,"data") 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]) ) + times <- seq(length(xdata))*stats::deltat(xdata) + stats::start(xdata)[1] + reconstructed <- lapply(seq(nfreq), function(i) ts( M$Amp[i]*cos(M$Freq[i]*times + M$Phase[i]), start=stats::start(xdata), deltat=stats::deltat(xdata)) ) + + if ( sum ) reconstructed <- ts(apply(simplify2array(reconstructed), 1 , sum), start=stats::start(xdata), deltat=stats::deltat(xdata)) } #' MFFT ANOVA diff --git a/R/toneCombinations.R b/R/toneCombinations.R new file mode 100644 index 0000000000000000000000000000000000000000..afcdab14d5c522c9d32f28532745e5e631eecf90 --- /dev/null +++ b/R/toneCombinations.R @@ -0,0 +1,113 @@ +#' Generation of combination of tones +#' +#' Generates a vector with combinations of an input vector of frequencies, wih +#' explicit label names, up to order 3 (this could be made more flexible is the future) +#' +#' @importFrom RcppAlgos +#' @param omegas: vector of references frequencies, optionally with rownames, +#' @param keepPositives : if TRUE, then only keeps positive combinations of frequencies +#' @return a vector with combination of tones and explicit rownames, using, if available, the +#' rownames provided in the input vector omega +#' @author Michel Crucifix +#' @examples +# omegas <- c( 0.123, 0.14312, 0.33251, 0.554313) +# outamps <- c(1., 2, 0.2 , 0.5, 0.5) +# outfreqs <- c(1., 1.2432, omegas[1]+omegas[3]+0.00000002, omegas[1]-omegas[4]+0.00004, 0.15) +# +# attributions <- attributeTones(outfreqs, omegas) +# +# cbind(outfreqs, attributions) +# +# plot(outfreqs, outamps, type='h') +# text(outfreqs, outamps+0.1, attributions) +# +toneCombinations <- function(omegas, keepPositives=TRUE){ + twoomegas <- c(-omegas,omegas) + indices <- c(-seq(length(omegas)), seq(length(omegas))) + result = rbind( + c(0, 0, 0), + cbind ( RcppAlgos::comboGeneral(twoomegas, m = 1, repetition = T), 0, 0), + cbind ( RcppAlgos::comboGeneral(twoomegas, m = 2, repetition = T) , 0), + RcppAlgos::comboGeneral(twoomegas, m = 3, repetition = T)) + combos = rbind( + c(0, 0, 0), + cbind ( RcppAlgos::comboGeneral(indices, m = 1, repetition = T), 0, 0), + cbind ( RcppAlgos::comboGeneral(indices, m = 2, repetition = T) , 0), + RcppAlgos::comboGeneral(indices, m = 3, repetition = T)) + + + sumre <- apply(result, 1, sum) + whichUnique <- seq(length(sumre))[!duplicated(sumre)] + uniqueCombos <- combos[whichUnique, ] + uniqueSums <- sumre[whichUnique ] + + # further filtering (?) + if (keepPositives){ + tokeep <- which(uniqueSums > 0) + uniqueCombos <- uniqueCombos[tokeep,] + uniqueSums <- uniqueSums[tokeep] + } + + names(uniqueSums) <- apply(uniqueCombos, 1, function(g) generate_name(g,"s", labels=names(omegas))) + return(cbind(uniqueSums)) + +} + +generate_name <- function(invec,char="s", labels = NULL){ + # to do: if labels are supplied, to not produce + tt <- as.data.frame(table(invec)) + flag <- FALSE + tmp <- sapply(seq(nrow(tt)), function(i) { + indice <- as.integer(as.character(tt[i,][[1]])) + if (indice == 0) return("") else { + if (flag) signchar <- ifelse(sign(indice)+1, " + "," - ") else + signchar <- ifelse(sign(indice)+1, "","-") + freq <- tt[i,][[2]] + freqchar <-ifelse(freq == 1, "", sprintf("%d",freq)) + flag <<- TRUE + if (is.null(labels)) labelname <- sprintf("%s%d", char, abs(indice)) else + labelname <- labels[abs(indice)] + return(sprintf("%s%s%s",signchar,freqchar,labelname)) + }}) +# Reduce(function(i) paste (i, sep=""), tmp) + Reduce(function(i, j) paste (i, j, sep=""), tmp) +} + + +#' Attribution of combination of tones +#' +#' Based on a vector of frequencies (`infreq`), and a vector of referenc +#' frequencies with row names (it will be input to `toneCombinations`), +#' attribute the `infreq` frequencies with two possible degrees of tolerance +#' +#' @param infreq : input frequencies +#' @param omegas : reference frequencies (a numeric vector which may contain explicit row names) +#' @param tol1 : acceptable tolerance for being considered as a certain attribution +#' (if several frequencies match the criteria, the closest will be taken) +#' @param tol2 : acceptable tolerance for being considered as a likely or plausible +#' +#' @examples +# omegas <- c( 0.123, 0.14312, 0.33251, 0.554313) +# names(omegas) <- c('g1','g2','s1','s2') +# outamps <- c(1., 2, 0.2 , 0.5, 0.5) +# outfreqs <- c(1., 1.2432, omegas[1]+omegas[3]+0.00000002, omegas[1]-omegas[4]+0.00004, 0.15) +# +# attributions <- attributeTones(outfreqs, omegas) +# +# cbind(outfreqs, attributions) +# +# plot(outfreqs, outamps, type='h') +# text(outfreqs, outamps+0.1, attributions) +# +attributeTones <- function(infreq , omegas, tol1 = 1.e-6, tol2 = 1.e-4) { + attributions <- rep("", length(infreq)) + combis <- toneCombinations(omegas) + for (i in seq(infreq)){ + deltas <- abs(infreq[i] - combis) + bestSuspect <- which.min(abs(infreq[i] - combis)) + if ( deltas[bestSuspect] < tol1 ) attributions[i] <- rownames(combis)[bestSuspect] else + if ( deltas[bestSuspect] < tol2 ) attributions[i] <- paste(rownames(combis)[bestSuspect],"?") + } + return (attributions) +} + diff --git a/man/mfft_deco.Rd b/man/mfft_deco.Rd index 0e6f5c9a8c92bfe999f502d085460edcd5a1cff2..1dd89634da5c781c9812f0c3fc97054b0a6c655c 100644 --- a/man/mfft_deco.Rd +++ b/man/mfft_deco.Rd @@ -9,7 +9,7 @@ \alias{print.mfft_deco} \title{MFFT reconstruction} \usage{ -reconstruct_mfft(M) +reconstruct_mfft(M, sum = TRUE) mfft_anova(M) @@ -21,6 +21,15 @@ mfft_anova(M) \method{print}{mfft_deco}(M, ...) } +\arguments{ +\item{M}{: mfft_deco object} + +\item{sum}{: TRUE if user wants to sum components in the reconstruction} +} +\value{ +list of reconstructed components if sum=FALSE, full + reconstructed time series otherwise +} \description{ MFFT reconstruction