From 3ce7a859271a373ae80ce05107dc872b461dc443 Mon Sep 17 00:00:00 2001
From: Michel Crucifix <michel.crucifix@uclouvain.be>
Date: Mon, 16 Dec 2024 09:01:28 +0100
Subject: [PATCH] fiddle with develop
---
R/develop.R | 32 ++++++++++++++++++++++++--------
R/mfft_complex.R | 9 +++++++--
R/mfft_real.R | 2 +-
3 files changed, 32 insertions(+), 11 deletions(-)
diff --git a/R/develop.R b/R/develop.R
index 3525a4e..0667212 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 78a3305..bd8e9d3 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 2f9fc87..488643b 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
--
GitLab