diff --git a/R/develop.R b/R/develop.R index 3525a4e3d39bdee00a86cc00ee4e95da3c31c016..0667212d7cb54ad19e0162fc5ccb2c991e29f129 100644 --- a/R/develop.R +++ b/R/develop.R @@ -32,15 +32,21 @@ cis <- function(x) exp(1i*x) #' reconstructed time series otherwise #' @method develop discreteSpectrum #' @export -develop.discreteSpectrum <- function(M, start=NULL, end=NULL, deltat=NULL, times, dfunction = cos, sum=TRUE){ +develop.discreteSpectrum <- function(M, start = NULL, end = NULL, deltat = NULL, times = NULL, dfunction = cos, maxfreq = NULL, sum=TRUE){ if (!("discreteSpectrum" %in% class(M))) stop ("object is not a discreteSpectrum decomposition") timesIsATseries = FALSE + if (is.ts(times)){ + start = start(times) + deltat= deltat(times) + end = start + length(times)*deltat + timesIsATseries <- TRUE + } if (!is.null(start)){ if (is.null(deltat) || is.null(end)) stop ("if you supply start, you must also supply deltat and end"); - n <- (end-start) %*% deltat - times <- start + seq(0, n) * deltat - timesIsATseries = TRUE + n <- (end-start) %/% deltat + times <- ts(start + seq(0, n) * deltat, start=start, deltat=deltat) + timesIsATseries <- TRUE } if (is.null(times)){ @@ -49,19 +55,29 @@ develop.discreteSpectrum <- function(M, start=NULL, end=NULL, deltat=NULL, time start <- stats::start(xdata)[1] deltat <- stats::deltat(xdata) times <- (seq(length(xdata))-1) * deltat + start - timesIsATseries = TRUE + timesIsATseries <- TRUE } nfreq <- attr(M,"nfreq") if (is.null(nfreq)) nfreq <- length(M$Amp) + if (!is.null(maxfreq)) nfreq <- min(nfreq, maxfreq) + if (timesIsATseries){ reconstructed <- lapply(seq(nfreq), function(i) ts( M$Amp[i] * dfunction(M$Freq[i] * times + M$Phase[i]), start=start, deltat=deltat) )} else { - reconstructed <- lapply(seq(nfreq), function(i) M$Amp[i] * dfunction(M$Freq[i] * times + M$Phase[i])) + reconstructed <- sapply(seq(nfreq), function(i) M$Amp[i] * dfunction(M$Freq[i] * times + M$Phase[i]) ) } - if ( sum ) reconstructed <- apply(simplify2array(reconstructed), 1 , sum) - if (timesIsATseries) reconstructed <- ts(reconstructed, start=start, deltat=deltat) +if ( sum ) { + shift <- attr(M, "shift"); if (is.null (shift)) shift = 0 + trend <- attr(M, "trend"); if (is.null (trend)) trend = 0 + if (timesIsATseries) + reconstructed <- Reduce('+', reconstructed) + trend * times + shift + else + reconstructed <- apply(reconstructed, 1 , sum) + trend * times + shift + } else if (timesIsATseries) { + reconstructed <- apply(reconstructed, 2, function(x) ts(x,start=start, deltat=deltat)) +} return(reconstructed) } diff --git a/R/mfft_complex.R b/R/mfft_complex.R index 78a3305bf4f39d8631fe8de7a0945aaa0a9e64ca..bd8e9d390df71cb9d927320ddef1fb32aa0fa60a 100644 --- a/R/mfft_complex.R +++ b/R/mfft_complex.R @@ -72,8 +72,13 @@ mfft_complex <- function(xdata, nfreq=30, minfreq=NULL, maxfreq=NULL, correctio xdata <- Re(xdata) ndata <- length(xdata); - flag <- c(1,'NA',2,3)[correction+1] - if (is.na(flag)) stop('this correction scheme is not implemented for complex time series') + flag <- c(1,0, 2,3)[correction+1] + print ('flag') + print (flag) + if (flag==0) { + message ('this correction scheme is not implemented for complex time series \\ will use synthetic data correction instead') + flag = 2 + } if (is.null(minfreq)){ my_minfreq <- -pi diff --git a/R/mfft_real.R b/R/mfft_real.R index 2f9fc873d5f7bce45e9b5049b263bf54133d408d..488643bbb0a0e1fe31e12fa46d299317eb4c4ce2 100644 --- a/R/mfft_real.R +++ b/R/mfft_real.R @@ -272,7 +272,7 @@ mfft_analyse <- function(xdata, nfreq, fast = TRUE, nu = NULL, minfreq=NULL, max # 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