-#Entrée : courbes synchrones, soit après étape 1 itérée, soit après chaqure étape 1
-#TODO: bout de code qui calcule les courbes synchrones après étapes 1+2 à partir des ID médoïdes
-
-#(Benjamin)
-#à partir de là, "conso" == courbes synchrones
-n <- nrow(conso)
-delta <- ncol(conso)
-
-#17000 colonnes coeff 1, puis 17000 coeff 2... [non : dans chaque tranche du cube]
-# #NOTE: delta et lscvect pourraient etre gardés à part (communs)
-
-#lignes 59 à 91 "dépliées" :
-Xcwt4 <- toCWT(conso, noctave = noctave4, dt = 1,
- scalevector = scalevector4,
- lt = delta, smooth = FALSE,
- nvoice = nvoice) # observations node with CWT
-
-#toCWT: (aux)
-##NOTE: renvoie une matrice 3D
- toCWT <- function(X, sw= 0, tw= 0, swabs= 0,
- nvoice= 12, noctave= 5,
- s0= 2, w0= 2*pi, lt= 24, dt= 0.5,
- spectra = FALSE, smooth = TRUE,
- scaled = FALSE,
- scalevector)
- { noctave <- adjust.noctave(lt, dt, s0, tw, noctave)
- if(missing(scalevector))
- scalevector <- 2^(0:(noctave * nvoice) / nvoice) * s0
- res <- lapply(1:nrow(X), function(n)
- { tsX <- ts( X[n,] )
- tsCent <- tsX - mean(tsX)
- if(scaled) tsCent <- ts(scale(tsCent))
- tsCent.cwt <- cwt.ts(tsCent, s0, noctave, nvoice, w0)
- tsCent.cwt
- } )
- if( spectra ) res <- lapply(res, function(l) Mod(l)^2 )
- if( smooth ) res <- lapply(res, smCWT, swabs = swabs,
- tw = tw, dt = dt,
- scalevector = scalevector)
- resArray <- array(NA, c(nrow(res[[1]]), ncol(res[[1]]),
- length(res)))
- for( l in 1:length(res) ) resArray[ , , l] <- res[[l]]
- resArray
- }
-
-#from sowas
-cwt.ts <- function(ts,s0,noctave=5,nvoice=10,w0=2*pi){
-
- if (class(ts)!="ts"){
-
- cat("# This function needs a time series object as input. You may construct this by using the function ts(data,start,deltat). Try '?ts' for help.\n")
-
- }
- else{
-
- t=time(ts)
- dt=t[2]-t[1]
-
- s0unit=s0/dt*w0/(2*pi)
- s0log=as.integer((log2(s0unit)-1)*nvoice+1.5)
-
- if (s0log<1){
- cat(paste("# s0unit = ",s0unit,"\n",sep=""))
- cat(paste("# s0log = ",s0log,"\n",sep=""))
- cat("# s0 too small for w0! \n")
- }
- totnoct=noctave+as.integer(s0log/nvoice)+1
-
- #cwt from package Rwave
- totts.cwt=cwt(ts,totnoct,nvoice,w0,plot=0)
-
- ts.cwt=totts.cwt[,s0log:(s0log+noctave*nvoice)]
-
- #Normalization
- sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0)
- smat <- matrix(rep(sqs,length(t)),nrow=length(t),byrow=TRUE)
-
- ts.cwt*smat
-
- }
-
+#precondition: ( log2(s0*w0/(2*pi)) - 1 ) * nvoice + 1.5 >= 1
+toCWT <- function(X, tw=0, swabs=0, nvoice=12, noctave=5, s0=2, w0=2*pi,
+ spectra=FALSE, smooth=TRUE, scaled=FALSE, scalevector)
+{
+ if(missing(scalevector))
+ scalevector <- 2^(0:(noctave * nvoice) / nvoice) * s0
+ s0log=as.integer((log2( s0*w0/(2*pi) )-1)*nvoice+1.5)
+ totnoct=noctave+as.integer(s0log/nvoice)+1
+ res <- lapply(1:nrow(X), function(n) {
+ ts <- scale(ts( X[n,] ), center=TRUE, scale=scaled)
+ totts.cwt = Rwave::cwt(ts,totnoct,nvoice,w0,plot=0)
+ ts.cwt=totts.cwt[,s0log:(s0log+noctave*nvoice)]
+ #Normalization
+ sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0)
+ smat <- matrix(rep(sqs,length(t)),nrow=length(t),byrow=TRUE)
+ ts.cwt*smat
+ })
+ if( spectra )
+ res <- lapply(res, function(l) Mod(l)^2 )
+ if( smooth )
+ res <- lapply(res, smCWT, swabs = swabs, tw = tw, scalevector = scalevector)
+ resArray <- array(NA, c(nrow(res[[1]]), ncol(res[[1]]), length(res)))
+ for( l in 1:length(res) )
+ resArray[ , , l] <- res[[l]]
+ resArray