| | 1 | #Ligne 83 dans epclust/R/computeWerDists.R : |
| | 2 | # TODO: |
| | 3 | # [X,period,scale,coix] = |
| | 4 | #wavelet(x(:,2),dt,Args.Pad,Args.Dj,Args.S0,Args.J1,Args.Mother);%#ok |
| | 5 | #[Y,period,scale,coiy] = wavelet(y(:,2),dt,Args.Pad,Args.Dj,Args.S0,Args.J1,Args.Mother); |
| | 6 | # |
| | 7 | #%Smooth X and Y before truncating! (minimize coi) |
| | 8 | #%sinv=1./(scale'); |
| | 9 | #% |
| | 10 | #% |
| | 11 | #%sX=smoothwavelet(sinv(:,ones(1,nx)).*(abs(X).^2),dt,period,Args.Dj,scale); |
| | 12 | #%sY=smoothwavelet(sinv(:,ones(1,ny)).*(abs(Y).^2),dt,period,Args.Dj,scale); |
| | 13 | #% |
| | 14 | #%Wxy=X.*conj(Y); |
| | 15 | #% |
| | 16 | #%% ----------------------- Wavelet coherence --------------------------------- |
| | 17 | #%sWxy=smoothwavelet(sinv(:,ones(1,n)).*Wxy,dt,period,Args.Dj,scale); |
| | 18 | #%Rsq=abs(sWxy).^2./(sX.*sY); |
| | 19 | # |
| | 20 | |
| | 21 | |
| | 22 | #https://github.com/grinsted/wavelet-coherence/blob/master/wtc.m |
| | 23 | |
| | 24 | |
| | 25 | library(Rwave) # CWT |
| | 26 | library(cluster) # pam |
| | 27 | #library(flexclust) # kcca |
| | 28 | |
| | 29 | |
| | 30 | toCWT <- function(X, sw= 0, tw= 0, swabs= 0, |
| | 31 | nvoice= 12, noctave= 5, |
| | 32 | s0= 2, w0= 2*pi, lt= 24, dt= 0.5, |
| | 33 | spectra = FALSE, smooth = TRUE, |
| | 34 | scaled = FALSE, |
| | 35 | scalevector) |
| | 36 | { noctave <- adjust.noctave(lt, dt, s0, tw, noctave) |
| | 37 | if(missing(scalevector)) |
| | 38 | scalevector <- 2^(0:(noctave * nvoice) / nvoice) * s0 |
| | 39 | res <- lapply(1:nrow(X), function(n) |
| | 40 | { tsX <- ts( X[n,] ) |
| | 41 | tsCent <- tsX - mean(tsX) |
| | 42 | if(scaled) tsCent <- ts(scale(tsCent)) |
| | 43 | tsCent.cwt <- cwt.ts(tsCent, s0, noctave, nvoice, w0) |
| | 44 | tsCent.cwt |
| | 45 | } ) |
| | 46 | if( spectra ) res <- lapply(res, function(l) Mod(l)^2 ) |
| | 47 | if( smooth ) res <- lapply(res, smCWT, swabs = swabs, |
| | 48 | tw = tw, dt = dt, |
| | 49 | scalevector = scalevector) |
| | 50 | resArray <- array(NA, c(nrow(res[[1]]), ncol(res[[1]]), |
| | 51 | length(res))) |
| | 52 | for( l in 1:length(res) ) resArray[ , , l] <- res[[l]] |
| | 53 | resArray |
| | 54 | } |
| | 55 | |
| | 56 | |
| | 57 | # =============================================================== |
| | 58 | |
| | 59 | smCWT <- function(CWT, sw= 0, tw= 0, swabs= 0, |
| | 60 | nvoice= 12, noctave= 2, s0= 2, w0= 2*pi, |
| | 61 | lt= 24, dt= 0.5, scalevector ) |
| | 62 | { |
| | 63 | # noctave <- adjust.noctave(lt, dt, s0, tw, noctave) |
| | 64 | # scalevector <- 2^(0:(noctave * nvoice) / nvoice) * s0 |
| | 65 | wsp <- Mod(CWT) |
| | 66 | smwsp <- smooth.matrix(wsp, swabs) |
| | 67 | smsmwsp <- smooth.time(smwsp, tw, dt, scalevector) |
| | 68 | smsmwsp |
| | 69 | } |
| | 70 | |
| | 71 | |
| | 72 | # =============================================================== |
| | 73 | |
| | 74 | toDWT <- function(x, filter.number = 6, family = "DaubLeAsymm") |
| | 75 | { x2 <- spline(x, n = 2^ceiling( log(length(x), 2) ), |
| | 76 | method = 'natural')$y |
| | 77 | Dx2 <- wd(x2, family = family, filter.number = filter.number)$D |
| | 78 | Dx2 |
| | 79 | } |
| | 80 | |
| | 81 | # =============================================================== |
| | 82 | |
| | 83 | contrib <- function(x) |
| | 84 | { J <- log( length(x)+1, 2) |
| | 85 | nrj <- numeric(J) |
| | 86 | t0 <- 1 |
| | 87 | t1 <- 0 |
| | 88 | for( j in 1:J ) { |
| | 89 | t1 <- t1 + 2^(J-j) |
| | 90 | nrj[j] <- sqrt( sum( x[t0:t1]^2 ) ) |
| | 91 | t0 <- t1 + 1 |
| | 92 | } |
| | 93 | return(nrj) |
| | 94 | } |
| | 95 | |
| | 96 | |
| | 97 | # ========================================= distance for coh === |
| | 98 | |
| | 99 | coherence <- function( x, y) |
| | 100 | { J <- log(length(x) + 1, 2) |
| | 101 | t0 <- 1 |
| | 102 | sg2_x <- 0 |
| | 103 | sg2_y <- 0 |
| | 104 | sg_xy <- 0 |
| | 105 | for(j in 0:(J - 1)) |
| | 106 | { t1 <- t0 + 2^(J - j)/2 - 1 |
| | 107 | tt <- t0:t1 |
| | 108 | sg2_x <- sg2_x + mean(x[t0:t1]^2) |
| | 109 | sg2_y <- sg2_y + mean(y[t0:t1]^2) |
| | 110 | sg_xy <- sg_xy + mean(x[t0:t1] * y[t0:t1]) |
| | 111 | t0 <- t1 + 1 |
| | 112 | } |
| | 113 | res <- sg_xy^2 / sg2_x / sg2_y |
| | 114 | res |
| | 115 | } |
| | 116 | |
| | 117 | |
| | 118 | vect2mat <- function(vect){ |
| | 119 | vect <- as.vector(vect) |
| | 120 | matrix(vect[-(1:2)], delta, lscvect) |
| | 121 | } |
| | 122 | |
| | 123 | |
| | 124 | # ========================================= # myimg for graphics |
| | 125 | jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", |
| | 126 | "cyan", "#7FFF7F", "yellow", |
| | 127 | "#FF7F00", "red", "#7F0000")) |
| | 128 | |
| | 129 | myimg <- function(MAT, x = 1:nrow(MAT), y = 1:col(MAT), ... ) |
| | 130 | filled.contour( x = x, y = y, z = MAT, |
| | 131 | xlab= 'Time', ylab= 'scale', |
| | 132 | color.palette = jet.colors, |
| | 133 | ... ) |
| | 134 | |
| | 135 | |
| | 136 | #TODO: [plus tard] alternative à sowa (package disparu) : cwt.. |
| | 137 | #source("sowas-superseded.r") # auxiliary CWT functions |
| | 138 | |
| | 139 | |
| | 140 | cwt.ts <- function(ts,s0,noctave=5,nvoice=10,w0=2*pi){ |
| | 141 | |
| | 142 | if (class(ts)!="ts"){ |
| | 143 | |
| | 144 | 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") |
| | 145 | |
| | 146 | } |
| | 147 | else{ |
| | 148 | |
| | 149 | t=time(ts) |
| | 150 | dt=t[2]-t[1] |
| | 151 | |
| | 152 | s0unit=s0/dt*w0/(2*pi) |
| | 153 | s0log=as.integer((log2(s0unit)-1)*nvoice+1.5) |
| | 154 | |
| | 155 | if (s0log<1){ |
| | 156 | cat(paste("# s0unit = ",s0unit,"\n",sep="")) |
| | 157 | cat(paste("# s0log = ",s0log,"\n",sep="")) |
| | 158 | cat("# s0 too small for w0! \n") |
| | 159 | } |
| | 160 | totnoct=noctave+as.integer(s0log/nvoice)+1 |
| | 161 | |
| | 162 | totts.cwt=cwt(ts,totnoct,nvoice,w0,plot=0) |
| | 163 | |
| | 164 | ts.cwt=totts.cwt[,s0log:(s0log+noctave*nvoice)] |
| | 165 | |
| | 166 | #Normalization |
| | 167 | sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0) |
| | 168 | smat <- matrix(rep(sqs,length(t)),nrow=length(t),byrow=TRUE) |
| | 169 | |
| | 170 | ts.cwt*smat |
| | 171 | |
| | 172 | } |
| | 173 | |
| | 174 | } |
| | 175 | |
| | 176 | ##################################### |
| | 177 | # WSP |
| | 178 | ##################################### |
| | 179 | |
| | 180 | wsp <- function(ts,s0=1,noctave=5,nvoice=10,w0=2*pi,sw=0,tw=0,swabs=0,siglevel=0.95,critval=NULL,nreal=1000,arealsiglevel=0.9,kernel=0,markt=-999,marks=-999,logscale=FALSE,plot=TRUE,units="",device="screen",file="wsp",color=TRUE,pwidth=10,pheight=7,labsc=1,labtext="",sigplot=3){ |
| | 181 | |
| | 182 | if (class(ts)!="ts"){ |
| | 183 | |
| | 184 | 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") |
| | 185 | |
| | 186 | } |
| | 187 | else{ |
| | 188 | |
| | 189 | if (sw!=0 & swabs==0) |
| | 190 | swabs <- as.integer(sw*nvoice) |
| | 191 | if (swabs!=0 & sw==0) |
| | 192 | sw <- swabs/nvoice |
| | 193 | |
| | 194 | sllist <- checkarealsiglevel(sw,tw,w0,arealsiglevel,siglevel,0) |
| | 195 | arealsiglevel <- sllist$arealsiglevel |
| | 196 | siglevel <- sllist$siglevel |
| | 197 | |
| | 198 | at <- NULL |
| | 199 | |
| | 200 | t <- time(ts) |
| | 201 | dt <- deltat(ts) |
| | 202 | s0rem <- s0 |
| | 203 | s0 <- adjust.s0(s0,dt) |
| | 204 | dnoctave <- as.integer(log(s0/s0rem-0.000000000001)/log(2))+1 |
| | 205 | |
| | 206 | noctave <- adjust.noctave(length(ts),dt,s0,tw,noctave) |
| | 207 | scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0 |
| | 208 | tsanom <- ts-mean(ts) |
| | 209 | |
| | 210 | #WAVELET TRANSFORMATION |
| | 211 | ts.cwt <- cwt.ts(tsanom,s0,noctave,nvoice,w0) |
| | 212 | |
| | 213 | #SMOOTHING |
| | 214 | wsp <- smooth.matrix(Mod(ts.cwt)^2,swabs) |
| | 215 | smwsp <- smooth.time(wsp,tw,dt,scalevector) |
| | 216 | |
| | 217 | #POINTWISE SIGNIFICANCE TEST |
| | 218 | if (is.null(critval)==FALSE){ # is critval empty? |
| | 219 | if (dim(critval)[2]!=dim(smwsp)[2]){ # is critval of the wrong dimension? |
| | 220 | if (siglevel[1]!=0 & nreal!=0) critval <- |
| | 221 | criticalvaluesWSP(tsanom,s0,noctave,nvoice,w0,swabs,tw,siglevel,nreal) |
| | 222 | #critval is of wrong dimension and siglevel and nreal are given |
| | 223 | else { |
| | 224 | critval <- NULL # no test possible |
| | 225 | arealsiglevel <- 0 |
| | 226 | cat("# dimension of critval is wrong \n") |
| | 227 | cat("# areawise test only possible with pointwise test \n") |
| | 228 | } |
| | 229 | } |
| | 230 | } |
| | 231 | else{ # critval is empty, but nreal or siglevel is given |
| | 232 | if (siglevel[1]!=0 & nreal!=0) critval <- |
| | 233 | criticalvaluesWSP(tsanom,s0,noctave,nvoice,w0,swabs,tw,siglevel,nreal) |
| | 234 | else { |
| | 235 | critval <- NULL |
| | 236 | arealsiglevel <- 0 |
| | 237 | cat("# areawise test only possible with pointwise test \n") |
| | 238 | } |
| | 239 | } |
| | 240 | |
| | 241 | #AREAL SIGNIFICANCE TEST |
| | 242 | if (arealsiglevel!=0){ |
| | 243 | v <- critval[1,] |
| | 244 | v[v==0] <- 9999 |
| | 245 | cvmat <- matrix(rep(v,length(t)),nrow=length(t),byrow=TRUE) |
| | 246 | atest <- arealtest(smwsp/cvmat,dt,s0,noctave,nvoice,w0,swabs,tw,siglevel,arealsiglevel,kernel,0) |
| | 247 | at <- atest$at |
| | 248 | kernel <- atest$kernel |
| | 249 | } |
| | 250 | |
| | 251 | if (s0rem<s0){ |
| | 252 | smwsp <- addvalues(nvoice,dnoctave,smwsp,NA) |
| | 253 | critval <- addvalues(nvoice,dnoctave,critval,1) |
| | 254 | #at <- addvalues(nvoice,dnoctave,at,NA) |
| | 255 | noctave <- noctave+dnoctave |
| | 256 | s0 <- s0/2^dnoctave |
| | 257 | scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0 |
| | 258 | } |
| | 259 | |
| | 260 | #PARAMETERS |
| | 261 | wclist <- |
| | 262 | list(modulus=smwsp,phase=NULL,time=t,s0=s0,noctave=noctave,nvoice=nvoice,w0=w0,scales=scalevector,critval=critval,at=at,kernel=kernel) |
| | 263 | |
| | 264 | class(wclist) <- "wt" |
| | 265 | |
| | 266 | #PLOTTING |
| | 267 | if (plot) |
| | 268 | plot(wclist,markt,marks,NULL,NULL,logscale,FALSE,units,"Wavelet Power Spectrum",device,file,FALSE,color,pwidth,pheight,labsc,labtext,sigplot) |
| | 269 | |
| | 270 | wclist |
| | 271 | |
| | 272 | } |
| | 273 | |
| | 274 | } |
| | 275 | |
| | 276 | ##################################### |
| | 277 | # WCO |
| | 278 | ##################################### |
| | 279 | |
| | 280 | wco <- |
| | 281 | function(ts1,ts2,s0=1,noctave=5,nvoice=10,w0=2*pi,sw=0,tw=0,swabs=0,siglevel=0.95,arealsiglevel=0.9,kernel=0,markt=-999,marks=-999,sqr=FALSE,phase=TRUE,plot=TRUE,units="",device="screen",file="wcoh",split=FALSE,color=TRUE,pwidth=10,pheight=7,labsc=1,labtext="",sigplot=3){ |
| | 282 | |
| | 283 | if (class(ts1)!="ts"){ |
| | 284 | |
| | 285 | cat("# This function needs two time series objects as input. You may construct them by using the function ts(data,start,deltat). Try '?ts' for help.\n") |
| | 286 | |
| | 287 | } |
| | 288 | else{ |
| | 289 | |
| | 290 | if (sw!=0 & swabs==0) |
| | 291 | swabs <- as.integer(sw*nvoice) |
| | 292 | if (swabs!=0 & sw==0) |
| | 293 | sw <- swabs/nvoice |
| | 294 | |
| | 295 | sllist <- checkarealsiglevel(sw,tw,w0,arealsiglevel,siglevel,1) |
| | 296 | arealsiglevel <- sllist$arealsiglevel |
| | 297 | siglevel <- sllist$siglevel |
| | 298 | |
| | 299 | if (sw==0 & tw==0 & swabs==0) { |
| | 300 | cat("# coherence without averaging makes no sense! \n") |
| | 301 | siglevel <- 0 |
| | 302 | arealsiglevel <- 0 |
| | 303 | } |
| | 304 | |
| | 305 | if (phase==FALSE) phs <- NULL |
| | 306 | |
| | 307 | at <- NULL |
| | 308 | |
| | 309 | tsadjust <- adjust.timeseries(ts1,ts2) |
| | 310 | ts1 <- tsadjust$ts1 |
| | 311 | ts2 <- tsadjust$ts2 |
| | 312 | |
| | 313 | t <- time(ts1) |
| | 314 | dt <- deltat(ts1) |
| | 315 | |
| | 316 | s0rem <- s0 |
| | 317 | s0 <- adjust.s0(s0,dt) |
| | 318 | dnoctave <- as.integer(log(s0/s0rem-0.000000000001)/log(2))+1 |
| | 319 | |
| | 320 | noctave <- adjust.noctave(length(ts1),dt,s0,tw,noctave) |
| | 321 | |
| | 322 | scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0 |
| | 323 | |
| | 324 | ts1anom <- ts1-mean(ts1) |
| | 325 | ts2anom <- ts2-mean(ts2) |
| | 326 | |
| | 327 | ts1.cwt <- cwt.ts(ts1anom,s0,noctave,nvoice,w0) |
| | 328 | ts2.cwt <- cwt.ts(ts2anom,s0,noctave,nvoice,w0) |
| | 329 | |
| | 330 | cosp <- Re(ts1.cwt)*Re(ts2.cwt) + Im(ts1.cwt)*Im(ts2.cwt) |
| | 331 | quad <- Im(ts1.cwt)*Re(ts2.cwt) - Re(ts1.cwt)*Im(ts2.cwt) |
| | 332 | wsp1 <- Mod(ts1.cwt)^2 |
| | 333 | wsp2 <- Mod(ts2.cwt)^2 |
| | 334 | |
| | 335 | smcosp <- smooth.matrix(cosp,swabs) |
| | 336 | smquad <- smooth.matrix(quad,swabs) |
| | 337 | smwsp1 <- smooth.matrix(wsp1,swabs) |
| | 338 | smwsp2 <- smooth.matrix(wsp2,swabs) |
| | 339 | |
| | 340 | smsmcosp <- smooth.time(smcosp,tw,dt,scalevector) |
| | 341 | smsmquad <- smooth.time(smquad,tw,dt,scalevector) |
| | 342 | smsmwsp1 <- smooth.time(smwsp1,tw,dt,scalevector) |
| | 343 | smsmwsp2 <- smooth.time(smwsp2,tw,dt,scalevector) |
| | 344 | |
| | 345 | if (sqr==FALSE) |
| | 346 | wcoh <- sqrt((smsmcosp^2+smsmquad^2)/(smsmwsp1*smsmwsp2)) |
| | 347 | else |
| | 348 | wcoh <- (smsmcosp^2+smsmquad^2)/(smsmwsp1*smsmwsp2) |
| | 349 | |
| | 350 | if (phase) |
| | 351 | phs <- atan2(smsmquad,smsmcosp) |
| | 352 | else phs <- NULL |
| | 353 | |
| | 354 | #POINTWISE SIGNIFICANCE TEST |
| | 355 | if (siglevel[1]!=0) critval <- criticalvaluesWCO(s0,noctave,nvoice,w0,swabs,tw,siglevel) |
| | 356 | else critval <- NULL |
| | 357 | if (sqr==TRUE & is.null(critval)==FALSE) |
| | 358 | critval <- critval^2 |
| | 359 | |
| | 360 | #AREAWISE SIGNIFICANCE TEST |
| | 361 | if (arealsiglevel!=0){ |
| | 362 | atest <- arealtest(wcoh/critval,dt,s0,noctave,nvoice,w0,swabs,tw,siglevel,arealsiglevel,kernel,1) |
| | 363 | at <- atest$at |
| | 364 | kernel <- atest$kernel |
| | 365 | } |
| | 366 | |
| | 367 | wcoh[1,1] <- 0 |
| | 368 | wcoh[1,2] <- 1 |
| | 369 | |
| | 370 | if (phase){ |
| | 371 | phs[1,1] <- -pi |
| | 372 | phs[1,2] <- pi |
| | 373 | } |
| | 374 | |
| | 375 | if (s0rem<s0){ |
| | 376 | wcoh <- addvalues(nvoice,dnoctave,wcoh,NA) |
| | 377 | phs <- addvalues(nvoice,dnoctave,phs,NA) |
| | 378 | noctave <- noctave+dnoctave |
| | 379 | s0 <- s0/2^dnoctave |
| | 380 | scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0 |
| | 381 | } |
| | 382 | |
| | 383 | wclist <- list(modulus=wcoh,phase=phs,s0=s0,noctave=noctave,nvoice=nvoice,w0=w0,time=t,scales=scalevector,critval=critval,at=at,kernel=kernel) |
| | 384 | |
| | 385 | class(wclist) <- "wt" |
| | 386 | |
| | 387 | if (plot) plot(wclist,markt,marks,NULL,NULL,FALSE,phase,units,"Wavelet Coherence",device,file,split,color,pwidth,pheight,labsc,labtext,sigplot) |
| | 388 | |
| | 389 | wclist |
| | 390 | |
| | 391 | } |
| | 392 | |
| | 393 | } |
| | 394 | |
| | 395 | ##################################### |
| | 396 | # WCS |
| | 397 | ##################################### |
| | 398 | |
| | 399 | wcs <- function(ts1,ts2,s0=1,noctave=5,nvoice=10,w0=2*pi,sw=0,tw=0,swabs=0,markt=-999,marks=-999,logscale=FALSE,phase=TRUE,plot=TRUE,units="",device="screen",file="wcsp",split=FALSE,color=TRUE,pwidth=10,pheight=7,labsc=1,labtext=""){ |
| | 400 | |
| | 401 | if (class(ts1)!="ts"){ |
| | 402 | |
| | 403 | cat("# This function needs two time series objects as input. You may construct them by using the function ts(data,start,deltat). Try '?ts' for help. \n") |
| | 404 | |
| | 405 | } |
| | 406 | else{ |
| | 407 | |
| | 408 | if (sw!=0 & swabs==0) |
| | 409 | swabs <- as.integer(sw*nvoice) |
| | 410 | if (swabs!=0 & sw==0) |
| | 411 | sw <- swabs/nvoice |
| | 412 | |
| | 413 | tsadjust <- adjust.timeseries(ts1,ts2) |
| | 414 | ts1 <- tsadjust$ts1 |
| | 415 | ts2 <- tsadjust$ts2 |
| | 416 | |
| | 417 | t <- time(ts1) |
| | 418 | dt <- deltat(ts1) |
| | 419 | |
| | 420 | s0rem <- s0 |
| | 421 | s0 <- adjust.s0(s0,dt) |
| | 422 | dnoctave <- as.integer(log(s0/s0rem-0.000000000001)/log(2))+1 |
| | 423 | |
| | 424 | noctave <- adjust.noctave(length(ts1),dt,s0,tw,noctave) |
| | 425 | |
| | 426 | scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0 |
| | 427 | |
| | 428 | ts1anom <- ts1-mean(ts1) |
| | 429 | ts2anom <- ts2-mean(ts2) |
| | 430 | |
| | 431 | ts1.cwt <- cwt.ts(ts1anom,s0,noctave,nvoice,w0) |
| | 432 | ts2.cwt <- cwt.ts(ts2anom,s0,noctave,nvoice,w0) |
| | 433 | |
| | 434 | cosp <- Re(ts1.cwt)*Re(ts2.cwt) + Im(ts1.cwt)*Im(ts2.cwt) |
| | 435 | quad <- Im(ts1.cwt)*Re(ts2.cwt) - Re(ts1.cwt)*Im(ts2.cwt) |
| | 436 | |
| | 437 | smcosp <- smooth.matrix(cosp,swabs) |
| | 438 | smquad <- smooth.matrix(quad,swabs) |
| | 439 | |
| | 440 | smsmcosp <- smooth.time(smcosp,tw,dt,scalevector) |
| | 441 | smsmquad <- smooth.time(smquad,tw,dt,scalevector) |
| | 442 | |
| | 443 | wcsp <- smsmcosp^2+smsmquad^2 |
| | 444 | |
| | 445 | if (phase) |
| | 446 | phs <- atan2(smsmquad,smsmcosp) |
| | 447 | else phs <- NULL |
| | 448 | |
| | 449 | if (phase){ |
| | 450 | phs[1,1] <- -pi |
| | 451 | phs[1,2] <- pi |
| | 452 | } |
| | 453 | |
| | 454 | if (s0rem<s0){ |
| | 455 | wcsp <- addvalues(nvoice,dnoctave,wcoh,NA) |
| | 456 | phs <- addvalues(nvoice,dnoctave,phs,NA) |
| | 457 | noctave <- noctave+dnoctave |
| | 458 | s0 <- s0/2^dnoctave |
| | 459 | scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0 |
| | 460 | } |
| | 461 | |
| | 462 | wclist <- list(modulus=wcsp,phase=phs,s0=s0,noctave=noctave,nvoice=nvoice,w0=w0,time=t,scales=scalevector,critval=NULL,at=NULL,kernel=NULL) |
| | 463 | |
| | 464 | class(wclist) <- "wt" |
| | 465 | |
| | 466 | if (plot) plot(wclist,markt,marks,NULL,NULL,logscale,phase,units,"Wavelet Cross Spectrum",device,file,split,color,pwidth,pheight,labsc,labtext) |
| | 467 | |
| | 468 | wclist |
| | 469 | |
| | 470 | } |
| | 471 | |
| | 472 | } |
| | 473 | |
| | 474 | ########################################## |
| | 475 | # POINTWISE SIGNIFICANCE TEST |
| | 476 | ########################################## |
| | 477 | |
| | 478 | rawWSP <- function(ts,s0=1,noctave=5,nvoice=20,w0=2*pi,swabs=0,tw=0,scalevector){ |
| | 479 | |
| | 480 | tsanom <- ts-mean(ts) |
| | 481 | |
| | 482 | ts.cwt <- cwt.ts(tsanom,s0,noctave,nvoice,w0) |
| | 483 | |
| | 484 | wsp <- Mod(ts.cwt)^2 |
| | 485 | |
| | 486 | smwsp <- smooth.matrix(wsp,swabs) |
| | 487 | smsmwsp <- smooth.time(smwsp,tw,deltat(ts),scalevector) |
| | 488 | |
| | 489 | smsmwsp |
| | 490 | |
| | 491 | } |
| | 492 | |
| | 493 | criticalvaluesWSP <- function(ts,s0=1,noctave=5,nvoice=10,w0=2*pi,swabs=0,tw=0,siglevel=0.95,nreal=1000){ |
| | 494 | |
| | 495 | t=time(ts) |
| | 496 | dt=deltat(ts) |
| | 497 | |
| | 498 | s0 <- adjust.s0(s0,dt) |
| | 499 | noctave <- adjust.noctave(length(ts),dt,s0,tw,noctave) |
| | 500 | |
| | 501 | scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0 |
| | 502 | |
| | 503 | cat("# calculating critical values by means of MC-simulations \n") |
| | 504 | |
| | 505 | s0unit=s0/dt*w0/(2*pi) |
| | 506 | s0log=as.integer((log2(s0unit)-1)*nvoice+1.5) |
| | 507 | |
| | 508 | wsp0 <- rawWSP(ts,s0,noctave,nvoice,w0,swabs,tw,scalevector) |
| | 509 | |
| | 510 | S <- dim(wsp0)[2] |
| | 511 | |
| | 512 | n1 <- 1 + 2*as.integer( sqrt(2) * 2^((S-swabs+s0log-1)/nvoice+1) ) |
| | 513 | n2 <- max(scalevector)*tw*2/dt*1.1 |
| | 514 | n3 <- 2*tw*s0*2^noctave/dt+100 |
| | 515 | n <- max(n1,n2,n3) |
| | 516 | |
| | 517 | center <- (n-1)/2+1 |
| | 518 | cv <- matrix(rep(0,S*length(siglevel)),ncol=S) |
| | 519 | rmatrix <- matrix(0,nrow=nreal,ncol=S) |
| | 520 | |
| | 521 | # Fitting AR1-process (red noise) to data |
| | 522 | arts0 <- ar(ts,order.max=1) |
| | 523 | sdts0 <- sd(ts[1:length(ts)]) |
| | 524 | |
| | 525 | if (arts0$order==0){ |
| | 526 | se <- sqrt(arts0$var) |
| | 527 | arts0$ar <- 0.000001 |
| | 528 | } |
| | 529 | else |
| | 530 | se <- sqrt(sdts0*sdts0*(1-arts0$ar*arts0$ar)) |
| | 531 | |
| | 532 | tsMC <- ts(data=rep(0,n),start=t[1],frequency=1/dt) |
| | 533 | |
| | 534 | # MC Simulations |
| | 535 | for (i in 1:nreal){ |
| | 536 | |
| | 537 | tsMC[1:n] <- arima.sim(list(ar = arts0$ar), n+1, sd = se)[2:(n+1)] |
| | 538 | |
| | 539 | rmatrix[i,] <- rawWSP(tsMC,s0,noctave,nvoice,w0,swabs,tw,scalevector)[center,] |
| | 540 | |
| | 541 | } |
| | 542 | |
| | 543 | for (s in (1+swabs):(S-swabs)) rmatrix[,s] <- sort(rmatrix[,s]) |
| | 544 | |
| | 545 | for (i in 1:length(siglevel)){ |
| | 546 | sigindex <- as.integer(nreal*siglevel[i]) |
| | 547 | cvv <- rmatrix[sigindex,] |
| | 548 | cvv[is.na(cvv)] <- 0 |
| | 549 | cv[i,] <- cvv |
| | 550 | } |
| | 551 | |
| | 552 | cv |
| | 553 | |
| | 554 | } |
| | 555 | |
| | 556 | ########### |
| | 557 | |
| | 558 | criticalvaluesWCO <- function(s0,noctave,nvoice,w0,swabs,tw,siglevel=0.95){ |
| | 559 | |
| | 560 | cv=rep(0,length(siglevel)) |
| | 561 | |
| | 562 | for (i in 1:length(siglevel)){ |
| | 563 | |
| | 564 | if (siglevel[i]!=0.9 && siglevel[i]!=0.95 && siglevel[i]!=0.99) siglevel[i] <- 0.95 |
| | 565 | |
| | 566 | if (siglevel[i]==0.9){ |
| | 567 | cat("# significance level set to 0.90 \n") |
| | 568 | sw <- 1.0*swabs/nvoice |
| | 569 | cv[i] <- 0.00246872*w0^2*sw + 0.0302419*w0*sw^2 + 0.342056*sw^3 - |
| | 570 | 0.000425853*w0^2 - 0.101749*w0*sw - 0.703537*sw^2 + |
| | 571 | 0.00816219*w0 + 0.442342*sw + 0.970315 |
| | 572 | } |
| | 573 | |
| | 574 | if (siglevel[i]==0.95){ |
| | 575 | cat("# significance level set to 0.95 \n") |
| | 576 | sw <- swabs*100.0/3.0/nvoice |
| | 577 | cv[i] <- 0.0000823*w0^3 + 0.0000424*w0^2*sw + 0.0000113*w0*sw^2 + |
| | 578 | 0.0000154*sw^3 - 0.0023*w0^2 - 0.00219*w0*sw - 0.000751*sw^2 + |
| | 579 | 0.0205*w0 + 0.0127*sw + 0.95 |
| | 580 | } |
| | 581 | |
| | 582 | if (siglevel[i]==0.99){ |
| | 583 | cat("# significance level set to 0.99 \n") |
| | 584 | sw <- 1.0*swabs/nvoice |
| | 585 | cv[i] <- 0.00102304*w0^2*sw + 0.00745772*w0*sw^2 + 0.230035*sw^3 - |
| | 586 | 0.000361565*w0^2 - 0.0502861*w0*sw - 0.440777*sw^2 + |
| | 587 | 0.00694998*w0 + 0.29643*sw + 0.972964 |
| | 588 | } |
| | 589 | |
| | 590 | if (cv[i]>1) cv[i] <- 1 |
| | 591 | |
| | 592 | cat(paste("# significance testing, cv=",cv[i]," \n",sep="")) |
| | 593 | |
| | 594 | } |
| | 595 | |
| | 596 | cv |
| | 597 | |
| | 598 | } |
| | 599 | |
| | 600 | ############################# |
| | 601 | # AREAWISE SIGNIFICANCE TEST |
| | 602 | ############################# |
| | 603 | |
| | 604 | slide <- function(data,kernellist,s0,noctave,nvoice,cv){ |
| | 605 | |
| | 606 | # slides kernel over data set |
| | 607 | #---------------------------- |
| | 608 | # data: matrix containing data |
| | 609 | # kernellist: matrix containing kernel |
| | 610 | # s0: lowest scale |
| | 611 | # noctave: number of octaves |
| | 612 | # nvoice: number of voices per octave |
| | 613 | # cv: critical value, all values higher are set to one |
| | 614 | |
| | 615 | #Time: (rows) n,i the kernel is to be scaled in this direction |
| | 616 | #Scale: (cols) m,j |
| | 617 | |
| | 618 | data[data<cv] <- 0 |
| | 619 | |
| | 620 | kernel <- kernellist$bitmap |
| | 621 | |
| | 622 | js <- kernellist$is |
| | 623 | |
| | 624 | sm <- s0*2^noctave |
| | 625 | |
| | 626 | dbin <- tobin(data) |
| | 627 | kbin <- tobin(kernel) |
| | 628 | |
| | 629 | dn <- nrow(dbin) |
| | 630 | dm <- ncol(dbin) |
| | 631 | |
| | 632 | kn <- nrow(kbin) |
| | 633 | km <- ncol(kbin) |
| | 634 | |
| | 635 | mark <- matrix(rep(0,dn*dm),nrow=dn) |
| | 636 | |
| | 637 | for (j in 1:(dm-km+1)){ |
| | 638 | |
| | 639 | s <- s0*2^((j+js-1)/nvoice) |
| | 640 | kscn <- as.integer(kn*s/sm); |
| | 641 | if (kscn==0) kscn <- 1 |
| | 642 | |
| | 643 | ksc <- scaleKernel(kbin,kscn) |
| | 644 | kscm <- km |
| | 645 | |
| | 646 | for (i in 1:(dn-kscn+1)){ |
| | 647 | |
| | 648 | subbin <- dbin[i:(i+kscn-1),j:(j+kscm-1)] |
| | 649 | |
| | 650 | if (sum(ksc*subbin)==sum(ksc)) |
| | 651 | mark[i:(i+kscn-1),j:(j+kscm-1)] <- mark[i:(i+kscn-1),j:(j+kscm-1)]+ksc |
| | 652 | |
| | 653 | } |
| | 654 | |
| | 655 | } |
| | 656 | |
| | 657 | mark <- tobin(mark) |
| | 658 | |
| | 659 | mark |
| | 660 | |
| | 661 | } |
| | 662 | |
| | 663 | arealtest <- function(wt,dt=1,s0=1,noctave=5,nvoice=20,w0=2*pi,swabs=0,tw=0,siglevel,arealsiglevel=0.9,kernel=0,type=0){ |
| | 664 | |
| | 665 | slp <- slope(w0,swabs,tw,nvoice,siglevel,arealsiglevel,type) |
| | 666 | |
| | 667 | if (length(kernel)<2){ |
| | 668 | maxarea <- s0*2^noctave*slp/10*nvoice/dt |
| | 669 | cat(paste("# calculating kernel (maxarea=",maxarea,")\n",sep="")) |
| | 670 | cvkernel <- |
| | 671 | kernelRoot(s0,w0,a=maxarea,noctave,nvoice,swabs,tw,dt) |
| | 672 | |
| | 673 | cat("# building kernel bitmap \n") |
| | 674 | kernel <- |
| | 675 | kernelBitmap(cvkernel,s0,w0,noctave,nvoice,swabs,tw,dt) |
| | 676 | |
| | 677 | } |
| | 678 | |
| | 679 | cat("# performing arealtest \n") |
| | 680 | sl <- slide(wt,kernel,s0,noctave,nvoice,1) |
| | 681 | |
| | 682 | list(at=sl,kernel=kernel) |
| | 683 | |
| | 684 | } |
| | 685 | |
| | 686 | ############################# |
| | 687 | # PLOTTING |
| | 688 | ############################# |
| | 689 | |
| | 690 | plotat <- function(t,wt,at,sigplot){ |
| | 691 | |
| | 692 | if (length(at)>1){ |
| | 693 | linewidth <- 1 |
| | 694 | if (sigplot==3) |
| | 695 | linewidth <- 5 |
| | 696 | |
| | 697 | contour(x=t,z=at,levels=0.5,add=TRUE,drawlabels=FALSE,axes=FALSE,lwd=linewidth,col="black") |
| | 698 | } |
| | 699 | |
| | 700 | } |
| | 701 | |
| | 702 | |
| | 703 | plotcv <- function(t,wt,cv){ |
| | 704 | |
| | 705 | if (length(dim(cv))==0) |
| | 706 | |
| | 707 | contour(x=t,z=wt,levels=c(cv),drawlabels=FALSE,axes=FALSE,add=TRUE,col="black",lwd=1) |
| | 708 | |
| | 709 | else{ |
| | 710 | |
| | 711 | for(i in 1:nrow(cv)){ |
| | 712 | |
| | 713 | v <- cv[i,] |
| | 714 | v[v==0] <- 9999 |
| | 715 | m <- matrix(rep(v,length(t)),nrow=length(t),byrow=TRUE) |
| | 716 | |
| | 717 | contour(x=t,z=wt/m,levels=1,drawlabels=FALSE,axes=FALSE,add=TRUE,col="black",lwd=1) |
| | 718 | |
| | 719 | } |
| | 720 | |
| | 721 | } |
| | 722 | |
| | 723 | } |
| | 724 | |
| | 725 | plotcoi <- function(t,s0,noctave,w0){ |
| | 726 | |
| | 727 | tv <- as.vector(t) |
| | 728 | tvl <- tv[tv-tv[1]<(tv[length(tv)]-tv[1])/2] |
| | 729 | tvr <- tv[tv-tv[1]>=(tv[length(tv)]-tv[1])/2] |
| | 730 | |
| | 731 | lines(tvl,log2(((tvl-tv[1])*4*pi/((w0+sqrt(2+w0*w0))*sqrt(2)))/s0)/noctave,col="black") |
| | 732 | lines(tvr,log2(((tv[length(tv)]-tvr)*4*pi/((w0+sqrt(2+w0*w0))*sqrt(2)))/s0)/noctave,col="black") |
| | 733 | |
| | 734 | } |
| | 735 | |
| | 736 | plotmarks <- function(t,s0,noctave,markt,marks){ |
| | 737 | |
| | 738 | if (markt[1]!=-999){ |
| | 739 | |
| | 740 | for (i in 1:length(markt)){ |
| | 741 | lines(c(markt[i],markt[i]),c(0,1),lty="dotted") |
| | 742 | } |
| | 743 | |
| | 744 | } |
| | 745 | |
| | 746 | if (marks[1]!=-999){ |
| | 747 | |
| | 748 | for (i in 1:length(marks)){ |
| | 749 | lines(c(t[1],t[length(t)]),c(log2(marks[i]/s0)/noctave,log2(marks[i]/s0)/noctave),lty="dotted") |
| | 750 | } |
| | 751 | |
| | 752 | } |
| | 753 | |
| | 754 | } |
| | 755 | |
| | 756 | |
| | 757 | ##################### |
| | 758 | # PLOT.WT |
| | 759 | ##################### |
| | 760 | |
| | 761 | plot.wt <- function(wt,markt=-999,marks=-999,t1=NULL,t2=NULL,logscale=FALSE,phase=FALSE,units="",plottitle="",device="screen",file="wt",split=FALSE,color=TRUE,pwidth=10,pheight=5,labsc=1.5,labtext="",sigplot=3,xax=NULL,xlab=NULL,yax=NULL,ylab=NULL){ |
| | 762 | |
| | 763 | plotwt(wt$modulus,wt$phase,wt$time,wt$s0,wt$noctave,wt$w0,wt$critval,wt$at,markt,marks,t1,t2,logscale,phase,units,plottitle,device,file,split,color,pwidth,pheight,labsc,labtext,sigplot,xax,xlab,yax,ylab) |
| | 764 | |
| | 765 | } |
| | 766 | |
| | 767 | plotwt <- |
| | 768 | function(wt,phs,t,s0,noctave,w0,cv=NULL,at=NULL,markt=-999,marks=-999,t1=NULL,t2=NULL,logscale=FALSE,phase=FALSE,units="",plottitle="Wavelet Plot",device="screen",file="wavelet",split=FALSE,color=TRUE,pwidth=10,pheight=7,labsc=1,labtext="",sigplot=1,xax=NULL,xlab=NULL,yax=NULL,ylab=NULL){ |
| | 769 | |
| | 770 | if (is.null(phs)) phase <- FALSE |
| | 771 | |
| | 772 | mgpv <- c(3,1,0) |
| | 773 | if (labsc>1) mgpv[1] <- 3-(labsc-1.5) |
| | 774 | |
| | 775 | ncolors <- 64 |
| | 776 | colors <- wtcolors(ncolors) |
| | 777 | if (color==FALSE) colors <- gray((0:ncolors)/ncolors*0.7+0.3) |
| | 778 | |
| | 779 | rangev <- (0:(ncolors-1))/(ncolors-1) |
| | 780 | rangebar <- matrix(rangev,nrow=2,ncol=64,byrow=TRUE) |
| | 781 | |
| | 782 | s <- 2^(0:(noctave))*s0 |
| | 783 | sn <- (0:(noctave))/noctave |
| | 784 | |
| | 785 | deltat <- (t[length(t)]-t[1])/(length(t)-1) |
| | 786 | cut <- FALSE |
| | 787 | if ((!is.null(t1)) | (!is.null(t2))){ |
| | 788 | if (t1<t2 & t1>=t[1] & t2<=t[length(t)]){ |
| | 789 | |
| | 790 | cut <- TRUE |
| | 791 | |
| | 792 | i1 <- (t1-t[1])/deltat+1 |
| | 793 | i2 <- (t2-t[1])/deltat+1 |
| | 794 | |
| | 795 | t <- t[t>=t1 & t<=t2] |
| | 796 | |
| | 797 | wt <- wt[i1:i2,] |
| | 798 | if (phase) phs <- phs[i1:i2,] |
| | 799 | if (length(at)>1) at <- at[i1:i2,] |
| | 800 | |
| | 801 | } |
| | 802 | } |
| | 803 | |
| | 804 | #---------------- |
| | 805 | # Setting layout |
| | 806 | #---------------- |
| | 807 | |
| | 808 | if (device=="ps" && split==FALSE){ |
| | 809 | |
| | 810 | file <- paste(file,".eps",sep="") |
| | 811 | |
| | 812 | postscript(file,onefile=TRUE,horizontal=TRUE,paper="special",width=pwidth,height=pheight) |
| | 813 | cat(paste("# writing",file," \n")) |
| | 814 | |
| | 815 | } |
| | 816 | |
| | 817 | if (device=="ps" && split==TRUE){ |
| | 818 | |
| | 819 | file1 <- paste(file,".mod.eps",sep="") |
| | 820 | |
| | 821 | postscript(file1,onefile=TRUE,horizontal=TRUE,paper="special",width=pwidth,height=pheight) |
| | 822 | cat(paste("# writing",file1," \n")) |
| | 823 | |
| | 824 | } |
| | 825 | |
| | 826 | if (phase==TRUE && split==FALSE) nlo <- layout(matrix(c(1,2,3,4),2,2,byrow=TRUE),widths=c(4,1)) |
| | 827 | else nlo <- layout(matrix(c(1,2),ncol=2,byrow=TRUE),widths=c(4,1)) |
| | 828 | |
| | 829 | |
| | 830 | #----------------- |
| | 831 | # Plotting modulus |
| | 832 | #----------------- |
| | 833 | |
| | 834 | if (logscale){ |
| | 835 | if (units=="") |
| | 836 | image(x=t,z=log10(wt),col = colors,axes=FALSE,xlab="Time",ylab="Scale",frame.plot=TRUE,cex.lab=labsc,mgp=mgpv) |
| | 837 | else |
| | 838 | image(x=t,z=log10(wt),col = colors,axes=FALSE,xlab=paste("Time ","[",units,"]",sep=""),ylab=paste("Scale ","[",units,"]",sep=""),frame.plot=TRUE,cex.lab=labsc,mgp=mgpv) |
| | 839 | } |
| | 840 | else{ |
| | 841 | if (units=="") |
| | 842 | image(x=t,z=wt,col = colors,axes=FALSE,xlab="Time",ylab="Scale",frame.plot=TRUE,cex.lab=labsc,mgp=mgpv) |
| | 843 | else |
| | 844 | image(x=t,z=wt,col = colors,axes=FALSE,xlab=paste("Time ","[",units,"]",sep=""),ylab=paste("Scale ","[",units,"]",sep=""),frame.plot=TRUE,cex.lab=labsc,mgp=mgpv) |
| | 845 | } |
| | 846 | |
| | 847 | text(t[1+as.integer(length(t)*0.1)],0.8,labtext,cex=2) |
| | 848 | |
| | 849 | if (sigplot==1 | sigplot==3){ |
| | 850 | if (is.null(cv)==FALSE){ #Critical values |
| | 851 | if (cv[1]!=0 & cv[1]!=-1) plotcv(t,wt,cv) |
| | 852 | } |
| | 853 | } |
| | 854 | |
| | 855 | if (sigplot>1) plotat(t,wt,at,sigplot) |
| | 856 | |
| | 857 | if (!cut) plotcoi(t,s0,noctave,w0) #cone of influence |
| | 858 | plotmarks(t,s0,noctave,markt,marks) #additional guiding lines |
| | 859 | |
| | 860 | if (is.null(xax)) |
| | 861 | axis(side=1,at=axTicks(1),cex.axis=labsc) |
| | 862 | else |
| | 863 | if (is.null(xlab)) |
| | 864 | axis(side=1,xax,labels=as.character(xax),cex.axis=labsc) |
| | 865 | else |
| | 866 | axis(side=1,xax,labels=xlab,cex.axis=labsc) |
| | 867 | |
| | 868 | if (is.null(yax)) |
| | 869 | axis(side=2,sn,labels=as.character(s),cex.axis=labsc) |
| | 870 | else |
| | 871 | if (is.null(ylab)) |
| | 872 | axis(side=2,yax,labels=as.character(yax),cex.axis=labsc) |
| | 873 | else |
| | 874 | axis(side=2,yax,labels=ylab,cex.axis=labsc) |
| | 875 | |
| | 876 | title(main=plottitle,cex=labsc) |
| | 877 | |
| | 878 | image(z=rangebar,axes=FALSE,col=colors,frame.plot=TRUE,cex.lab=labsc,mgp=mgpv) |
| | 879 | |
| | 880 | if (is.null(cv)==FALSE){ |
| | 881 | if (length(dim(cv))==0){ |
| | 882 | for (i in 1:length(cv)) |
| | 883 | if (cv[i]!=0) lines(c(-1,2),c(cv[i],cv[i])) |
| | 884 | } |
| | 885 | } |
| | 886 | |
| | 887 | if (!logscale) |
| | 888 | axis(side=2,(0:5)/5,labels=c("0","","","","","1"),cex.axis=labsc) |
| | 889 | else{ |
| | 890 | labelv <- substr(as.character(c(0:5)*(max(log10(wt),na.rm=TRUE)-min(log10(wt),na.rm=TRUE))/5),1,4) |
| | 891 | axis(side=2,(0:5)/5,labels=labelv,cex.axis=labsc) |
| | 892 | } |
| | 893 | |
| | 894 | |
| | 895 | #----------------- |
| | 896 | # Plotting phase |
| | 897 | #----------------- |
| | 898 | if (phase==TRUE){ |
| | 899 | |
| | 900 | if (device=="ps" && split==TRUE){ |
| | 901 | |
| | 902 | dev.off() |
| | 903 | |
| | 904 | file2 <- paste(file,".phs.eps",sep="") |
| | 905 | |
| | 906 | postscript(file2,onefile=TRUE,horizontal=TRUE,paper="special",width=10,height=5) |
| | 907 | cat(paste("# writing",file2," \n")) |
| | 908 | |
| | 909 | } |
| | 910 | |
| | 911 | colors <- rainbow(ncolors) |
| | 912 | if (color==FALSE){ |
| | 913 | dummy <- gray((0:ncolors)/ncolors) |
| | 914 | colors[1:(ncolors/2)] <- dummy[(ncolors/2+1):ncolors] |
| | 915 | colors[(ncolors/2+1):ncolors] <- dummy[1:(ncolors/2)] |
| | 916 | } |
| | 917 | |
| | 918 | if (units=="") |
| | 919 | image(x=t,z=phs,col=colors,axes=FALSE,xlab="Time",ylab="Scale",frame.plot=TRUE,cex.lab=labsc,mgp=mgpv) |
| | 920 | else |
| | 921 | image(x=t,z=phs,col=colors,axes=FALSE,xlab=paste("Time ","[",units,"]",sep=""),ylab=paste("Scale ","[",units,"]",sep=""),frame.plot=TRUE,cex.lab=labsc,mgp=mgpv) |
| | 922 | |
| | 923 | if (is.null(cv)==FALSE) plotcv(t,wt,cv) |
| | 924 | if (sigplot>1) plotat(t,wt,at,sigplot) |
| | 925 | if (!cut) plotcoi(t,s0,noctave,w0) |
| | 926 | plotmarks(t,s0,noctave,markt,marks) |
| | 927 | |
| | 928 | if (is.null(xax)) |
| | 929 | axis(side=1,at=axTicks(1),cex.axis=labsc) |
| | 930 | else |
| | 931 | if (is.null(xlab)) |
| | 932 | axis(side=1,xax,labels=as.character(xax),cex.axis=labsc) |
| | 933 | else |
| | 934 | axis(side=1,xax,labels=xlab,cex.axis=labsc) |
| | 935 | |
| | 936 | if (is.null(yax)) |
| | 937 | axis(side=2,sn,labels=as.character(s),cex.axis=labsc) |
| | 938 | else |
| | 939 | if (is.null(ylab)) |
| | 940 | axis(side=2,yax,labels=as.character(yax),cex.axis=labsc) |
| | 941 | else |
| | 942 | axis(side=2,yax,labels=ylab,cex.axis=labsc) |
| | 943 | |
| | 944 | |
| | 945 | title(main="Phase") |
| | 946 | |
| | 947 | image(z=rangebar,axes=FALSE,col=colors,frame.plot=TRUE,cex.lab=labsc,mgp=mgpv) |
| | 948 | axis(side=2,(0:4)/4,labels=c("-PI","","0","","PI"),cex.axis=labsc) |
| | 949 | |
| | 950 | } |
| | 951 | |
| | 952 | if (device=="ps") dev.off() |
| | 953 | |
| | 954 | } |
| | 955 | |
| | 956 | ############################## |
| | 957 | # Surrogates |
| | 958 | ############################## |
| | 959 | |
| | 960 | createwavesurrogates <- function(nsurr=1,wt=1,n,dt=1,s0=1,noctave=5,nvoice=10,w0=2*pi){ |
| | 961 | |
| | 962 | surrmat <- matrix(rep(0,n*nsurr),ncol=nsurr) |
| | 963 | |
| | 964 | for (i in 1:nsurr){ |
| | 965 | |
| | 966 | cat(paste("# Creating realization ",i,"\n",sep="")) |
| | 967 | |
| | 968 | x <- rnorm(n) |
| | 969 | xts <- ts(x,start=0,deltat=dt) |
| | 970 | |
| | 971 | xwt <- cwt.ts(xts,s0,noctave,nvoice,w0) |
| | 972 | wtsurr <- wt*xwt |
| | 973 | |
| | 974 | surri <- as.vector(invmorlet(wtsurr,0,dt,s0,noctave,nvoice,w0)) |
| | 975 | |
| | 976 | surrmat[,i] <- Re(surri) |
| | 977 | |
| | 978 | } |
| | 979 | |
| | 980 | surrmat |
| | 981 | |
| | 982 | } |
| | 983 | |
| | 984 | surrwave <- function(x,...) |
| | 985 | UseMethod("surrwave") |
| | 986 | |
| | 987 | surrwave.wt <- function(wt,nsurr=1,spec=TRUE){ |
| | 988 | |
| | 989 | n <- length(wt$time) |
| | 990 | t0 <- wt$time[1] |
| | 991 | dt <- (wt$time[13]-t0)/12 |
| | 992 | s0 <- wt$s0 |
| | 993 | noctave <- wt$noctave |
| | 994 | nvoice <- wt$nvoice |
| | 995 | w0 <- wt$w0 |
| | 996 | |
| | 997 | wt <- wt$modulus |
| | 998 | if (spec==TRUE) |
| | 999 | wt <- sqrt(wt) |
| | 1000 | |
| | 1001 | surrmat <- createwavesurrogates(nsurr,wt,n,dt,s0,noctave,nvoice,w0) |
| | 1002 | |
| | 1003 | surrts <- ts(surrmat,start=t0,deltat=dt) |
| | 1004 | |
| | 1005 | surrts |
| | 1006 | |
| | 1007 | } |
| | 1008 | |
| | 1009 | surrwave.matrix <- function(mat,nsurr=1,t0=0,dt=1,s0=1,noctave=5,nvoice=10,w0=2*pi,sw=0,tw=0,swabs=0,spec=FALSE){ |
| | 1010 | |
| | 1011 | if (sw!=0 & swabs==0) |
| | 1012 | swabs <- as.integer(sw*nvoice) |
| | 1013 | |
| | 1014 | scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0 |
| | 1015 | |
| | 1016 | if ((noctave*nvoice+1)!=dim(mat)[2]) |
| | 1017 | cat("# ERROR! nscales unequal noctave*nvoice+1 ! \n") |
| | 1018 | |
| | 1019 | n <- dim(mat)[1] |
| | 1020 | |
| | 1021 | if (spec==FALSE) |
| | 1022 | mat <- Mod(mat) |
| | 1023 | else |
| | 1024 | mat <- sqrt(Mod(mat)) |
| | 1025 | |
| | 1026 | wtsm <- smooth.matrix(mat,swabs) |
| | 1027 | wt <- smooth.time(wtsm,tw,dt,scalevector) |
| | 1028 | |
| | 1029 | surrmat <- createwavesurrogates(nsurr,wt,n,dt,s0,noctave,nvoice,w0) |
| | 1030 | |
| | 1031 | surrts <- ts(surrmat,start=t0,deltat=dt) |
| | 1032 | |
| | 1033 | surrts |
| | 1034 | |
| | 1035 | } |
| | 1036 | |
| | 1037 | surrwave.character <- |
| | 1038 | function(file,nsurr=1,t0=0,dt=1,s0=1,noctave=5,nvoice=10,w0=2*pi,sw=0,tw=0,swabs=0,transpose=TRUE,spec=FALSE){ |
| | 1039 | |
| | 1040 | if (sw!=0 & swabs==0) |
| | 1041 | swabs <- as.integer(sw*nvoice) |
| | 1042 | |
| | 1043 | scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0 |
| | 1044 | |
| | 1045 | if (transpose==FALSE) |
| | 1046 | mat <- matrix(scan(file,comment.char="#"),ncol=nvoice*noctave+1,byrow=TRUE) |
| | 1047 | else |
| | 1048 | mat <- matrix(scan(file,comment.char="#"),ncol=nvoice*noctave+1,byrow=FALSE) |
| | 1049 | |
| | 1050 | if ((noctave*nvoice+1)!=dim(mat)[2]) |
| | 1051 | cat("# ERROR! nscales unequal noctave*nvoice+1 ! \n") |
| | 1052 | |
| | 1053 | n <- dim(mat)[1] |
| | 1054 | |
| | 1055 | if (spec==FALSE) |
| | 1056 | mat <- Mod(mat) |
| | 1057 | else |
| | 1058 | mat <- sqrt(Mod(mat)) |
| | 1059 | |
| | 1060 | wtsm <- smooth.matrix(mat,swabs) |
| | 1061 | wt <- smooth.time(wtsm,tw,dt,scalevector) |
| | 1062 | |
| | 1063 | surrmat <- createwavesurrogates(nsurr,wt,n,dt,s0,noctave,nvoice,w0) |
| | 1064 | |
| | 1065 | surrts <- ts(surrmat,start=t0,deltat=dt) |
| | 1066 | |
| | 1067 | surrts |
| | 1068 | |
| | 1069 | } |
| | 1070 | |
| | 1071 | surrwave.ts <- function(ts,nsurr=1,s0=1,noctave=5,nvoice=10,w0=2*pi,sw=0,tw=0,swabs=0){ |
| | 1072 | |
| | 1073 | n <- length(ts) |
| | 1074 | t0 <- time(ts)[1] |
| | 1075 | dt <- deltat(ts) |
| | 1076 | if (sw!=0 & swabs==0) |
| | 1077 | swabs <- as.integer(sw*nvoice) |
| | 1078 | |
| | 1079 | scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0 |
| | 1080 | |
| | 1081 | wt <- Mod(cwt.ts(ts,s0,noctave,nvoice,w0)) |
| | 1082 | |
| | 1083 | wtsm <- smooth.matrix(wt,swabs) |
| | 1084 | wt <- smooth.time(wtsm,tw,dt,scalevector) |
| | 1085 | |
| | 1086 | surrmat <- createwavesurrogates(nsurr,wt,n,dt,s0,noctave,nvoice,w0) |
| | 1087 | |
| | 1088 | surrts <- ts(surrmat,start=t0,deltat=dt) |
| | 1089 | |
| | 1090 | surrts |
| | 1091 | |
| | 1092 | } |
| | 1093 | |
| | 1094 | invmorlet <- function(wt,t0=0,dt=1,s0=1,noctave=5,nvoice=10,w0=2*pi){ |
| | 1095 | |
| | 1096 | if ((noctave*nvoice+1)!=dim(wt)[2]) |
| | 1097 | cat("# ERROR! nscales unequal noctave*nvoice+1 ! \n") |
| | 1098 | |
| | 1099 | n <- dim(wt)[1] |
| | 1100 | |
| | 1101 | wt[is.na(wt)] <- 0 |
| | 1102 | |
| | 1103 | tsRe <- rep(0,n) |
| | 1104 | tsIm <- rep(0,n) |
| | 1105 | |
| | 1106 | wtRe <- t(Re(wt)) |
| | 1107 | wtIm <- t(Im(wt)) |
| | 1108 | |
| | 1109 | z <- .C("invmorlet", |
| | 1110 | as.double(as.vector(wtRe)), |
| | 1111 | as.double(as.vector(wtIm)), |
| | 1112 | as.integer(n), |
| | 1113 | as.double(dt), |
| | 1114 | as.double(s0), |
| | 1115 | as.integer(noctave), |
| | 1116 | as.integer(nvoice), |
| | 1117 | as.double(w0), |
| | 1118 | tsRetmp = as.double(tsRe), |
| | 1119 | tsImtmp = as.double(tsIm), |
| | 1120 | PACKAGE="sowas") |
| | 1121 | |
| | 1122 | invvec=complex(real=z$tsRetmp,imaginary=z$tsImtmp) |
| | 1123 | invts <- ts(data=invvec,start=t0,deltat=dt) |
| | 1124 | |
| | 1125 | invts |
| | 1126 | |
| | 1127 | } |
| | 1128 | |
| | 1129 | ################################# |
| | 1130 | # INPUT / OUTPUT |
| | 1131 | ################################# |
| | 1132 | |
| | 1133 | readmatrix <- function(file,M){ |
| | 1134 | |
| | 1135 | A <- matrix(scan(file,comment.char="#"),ncol=M,byrow=TRUE) |
| | 1136 | |
| | 1137 | A |
| | 1138 | |
| | 1139 | } |
| | 1140 | |
| | 1141 | |
| | 1142 | readts <- function(file){ |
| | 1143 | |
| | 1144 | A <- matrix(scan(file,comment.char="#"),ncol=2,byrow=TRUE) |
| | 1145 | |
| | 1146 | Adum <- A |
| | 1147 | |
| | 1148 | Adum[is.na(Adum)] <- 0 |
| | 1149 | |
| | 1150 | t <- Adum %*% c(1,0) |
| | 1151 | x <- A %*% c(0,1) |
| | 1152 | |
| | 1153 | N=length(t) |
| | 1154 | f=1/(t[13]-t[1])*12 |
| | 1155 | |
| | 1156 | if ((f>11) && (f<13)) f <- 12 |
| | 1157 | |
| | 1158 | timeseries<-ts(data=x,start=t[1],frequency=f) |
| | 1159 | |
| | 1160 | timeseries |
| | 1161 | |
| | 1162 | } |
| | 1163 | |
| | 1164 | writematrix <- function(file,data,header="# R Matrix"){ |
| | 1165 | |
| | 1166 | write(header,file) |
| | 1167 | write(t(data),file,ncol(data),append=TRUE) |
| | 1168 | |
| | 1169 | } |
| | 1170 | |
| | 1171 | ############################ |
| | 1172 | # HELP FUNCTIONS |
| | 1173 | ############################ |
| | 1174 | |
| | 1175 | smooth.matrix <- function(wt,swabs){ |
| | 1176 | |
| | 1177 | if (swabs != 0) |
| | 1178 | smwt <- t(filter(t(wt),rep(1,2*swabs+1)/(2*swabs+1))) |
| | 1179 | else |
| | 1180 | smwt <- wt |
| | 1181 | |
| | 1182 | smwt |
| | 1183 | |
| | 1184 | } |
| | 1185 | |
| | 1186 | smooth.time <- function(wt,tw,dt,scalevector){ |
| | 1187 | |
| | 1188 | smwt <- wt |
| | 1189 | |
| | 1190 | if (tw != 0){ |
| | 1191 | for (i in 1:length(scalevector)){ |
| | 1192 | |
| | 1193 | twi <- as.integer(scalevector[i]*tw/dt) |
| | 1194 | smwt[,i] <- filter(wt[,i],rep(1,2*twi+1)/(2*twi+1)) |
| | 1195 | |
| | 1196 | } |
| | 1197 | } |
| | 1198 | |
| | 1199 | smwt |
| | 1200 | |
| | 1201 | } |
| | 1202 | |
| | 1203 | |
| | 1204 | adjust.noctave <- function(N,dt,s0,tw,noctave){ |
| | 1205 | |
| | 1206 | if (tw>0){ |
| | 1207 | dumno <- as.integer((log(N*dt)-log(2*tw*s0))/log(2)) |
| | 1208 | if (dumno<noctave){ |
| | 1209 | cat("# noctave adjusted to time smoothing window \n") |
| | 1210 | noctave <- dumno |
| | 1211 | } |
| | 1212 | } |
| | 1213 | |
| | 1214 | noctave |
| | 1215 | |
| | 1216 | } |
| | 1217 | |
| | 1218 | adjust.s0 <- function(s0,dt){ |
| | 1219 | |
| | 1220 | if (s0<2*dt){ |
| | 1221 | s0 <- 2*dt |
| | 1222 | cat(paste("# s0 set to ",s0," \n")) |
| | 1223 | } |
| | 1224 | |
| | 1225 | s0 |
| | 1226 | |
| | 1227 | } |
| | 1228 | |
| | 1229 | adjust.timeseries <- function(ts1,ts2){ |
| | 1230 | |
| | 1231 | if (length(ts1)!=length(ts2)){ |
| | 1232 | tsint <- ts.intersect(ts1,ts2) |
| | 1233 | dt <- deltat(ts1) |
| | 1234 | ts1 <- ts(data=tsint[,1],start=time(tsint)[1],frequency=1/dt) |
| | 1235 | ts2 <- ts(data=tsint[,2],start=time(tsint)[1],frequency=1/dt) |
| | 1236 | t <- time(ts1) |
| | 1237 | } |
| | 1238 | |
| | 1239 | list(ts1=ts1,ts2=ts2) |
| | 1240 | |
| | 1241 | } |
| | 1242 | |
| | 1243 | checkarealsiglevel <- function(sw,tw,w0,arealsiglevel,siglevel,type){ |
| | 1244 | |
| | 1245 | if (type==0){ |
| | 1246 | |
| | 1247 | swv <- c(0,0.5,1) |
| | 1248 | twv <- c(0,1.5,3) |
| | 1249 | w0v <- c(pi,2*pi) |
| | 1250 | |
| | 1251 | if (length(swv[swv==sw])==0 || length(twv[twv==tw])==0 || |
| | 1252 | length(w0v[w0v==w0])==0){ |
| | 1253 | arealsiglevel <- 0 |
| | 1254 | cat("# areawise test for spectrum currently \n") |
| | 1255 | cat("# only possible for \n") |
| | 1256 | cat("# sw = 0 \n") |
| | 1257 | cat("# tw = 0 \n") |
| | 1258 | cat("# w0 = 2pi \n") |
| | 1259 | cat("# No areawise test performed \n") |
| | 1260 | } |
| | 1261 | } |
| | 1262 | |
| | 1263 | if (type==1){ |
| | 1264 | |
| | 1265 | swv <- c(0.5) |
| | 1266 | twv <- c(1.5) |
| | 1267 | w0v <- c(2*pi) |
| | 1268 | |
| | 1269 | if (length(swv[swv==sw])==0 || length(twv[twv==tw])==0 || |
| | 1270 | length(w0v[w0v==w0])==0){ |
| | 1271 | arealsiglevel <- 0 |
| | 1272 | cat("# areawise test for coherence currently \n") |
| | 1273 | cat("# only possible for \n") |
| | 1274 | cat("# sw = 0.5 \n") |
| | 1275 | cat("# tw = 1.5 \n") |
| | 1276 | cat("# w0 = 2pi \n") |
| | 1277 | cat("# No areawise test performed \n") |
| | 1278 | } |
| | 1279 | } |
| | 1280 | |
| | 1281 | |
| | 1282 | if (arealsiglevel!=0){ |
| | 1283 | arealsiglevel <- 0.9 |
| | 1284 | siglevel <- 0.95 |
| | 1285 | cat("# currently only siglevel=0.95 and arealsiglevel=0.9 possible for areawise test \n") |
| | 1286 | } |
| | 1287 | |
| | 1288 | list(siglevel=siglevel,arealsiglevel=arealsiglevel) |
| | 1289 | |
| | 1290 | } |
| | 1291 | |
| | 1292 | ######################## |
| | 1293 | |
| | 1294 | as.wt <- function(modulus,phase=NULL,s0=NULL,noctave=NULL,nvoice=NULL,w0=NULL,dt=1,time=NULL,scales=NULL,critval=NULL,at=NULL,kernel=NULL,N=NULL,t0=NULL){ |
| | 1295 | |
| | 1296 | if (is.null(scales)) |
| | 1297 | gotscales <- FALSE |
| | 1298 | else |
| | 1299 | gotscales <- TRUE |
| | 1300 | |
| | 1301 | if ((!gotscales) & (!is.null(s0)) & (!is.null(noctave)) &(!is.null(nvoice))){ |
| | 1302 | gotscales <- TRUE |
| | 1303 | scales=2^(0:(noctave*nvoice)/nvoice)*s0 |
| | 1304 | } |
| | 1305 | |
| | 1306 | if (gotscales & (is.null(s0) | is.null(noctave) | |
| | 1307 | is.null(nvoice))){ |
| | 1308 | s0 <- scales[1] |
| | 1309 | noctave <- log(scales[length(scales)]/s0)/log(2) |
| | 1310 | nvoice <- (length(scales)-1)/noctave |
| | 1311 | } |
| | 1312 | |
| | 1313 | |
| | 1314 | if (!gotscales) |
| | 1315 | cat("# ERROR! No scales given! \n") |
| | 1316 | |
| | 1317 | if (is.null(time)) gottimes <- FALSE |
| | 1318 | else gottimes <- TRUE |
| | 1319 | |
| | 1320 | if ((!gottimes) & (!is.null(dt)) & (!is.null(t0)) &(!is.null(N))){ |
| | 1321 | gottimes <- TRUE |
| | 1322 | time=(0:(N-1))*dt+t0 |
| | 1323 | } |
| | 1324 | |
| | 1325 | if (!gottimes) |
| | 1326 | cat("# ERROR! No time vector given! \n") |
| | 1327 | |
| | 1328 | wcolist <- list(modulus=modulus,phase=phase,s0=s0,noctave=noctave,nvoice=nvoice,w0=w0,time=time,scales=scales,critval=critval,at=at,kernel=kernel) |
| | 1329 | |
| | 1330 | class(wcolist) <- "wt" |
| | 1331 | |
| | 1332 | wcolist |
| | 1333 | |
| | 1334 | } |
| | 1335 | |
| | 1336 | ######################## |
| | 1337 | |
| | 1338 | wtcolors <- function(ncolors){ |
| | 1339 | |
| | 1340 | upside <- rainbow(ncolors,start=0,end=.7) |
| | 1341 | #upside <- heat.colors(ncolors+5) |
| | 1342 | #upside <- upside[1:ncolors] |
| | 1343 | |
| | 1344 | |
| | 1345 | down <- upside |
| | 1346 | |
| | 1347 | for (i in 1:ncolors){ |
| | 1348 | down[i] <- upside[ncolors+1-i] |
| | 1349 | }# |
| | 1350 | |
| | 1351 | down |
| | 1352 | |
| | 1353 | } |
| | 1354 | |
| | 1355 | #################### |
| | 1356 | |
| | 1357 | createwgn <- function(N,sig,dt){ |
| | 1358 | |
| | 1359 | timeseries<-ts(data=rnorm(N,0,sig),start=0,deltat=dt) |
| | 1360 | |
| | 1361 | timeseries |
| | 1362 | |
| | 1363 | } |
| | 1364 | |
| | 1365 | |
| | 1366 | createar <- function(N,a,sig,dt){ |
| | 1367 | |
| | 1368 | if (a==0) a <- 0.000000001 |
| | 1369 | |
| | 1370 | se <- sqrt(sig*sig*(1-a*a)) |
| | 1371 | tsMC <- ts(data=rep(0,N),start=0,deltat=dt) |
| | 1372 | |
| | 1373 | tsMC[1:N] <- arima.sim(list(ar = a), N, sd = se) |
| | 1374 | |
| | 1375 | } |
| | 1376 | |
| | 1377 | ###################### |
| | 1378 | |
| | 1379 | rk <- function(N=1000,s=8,noctave=5,nvoice=10,w0=2*pi,plot=TRUE){ |
| | 1380 | |
| | 1381 | t <- 1:N |
| | 1382 | |
| | 1383 | sunit <- s*(w0+sqrt(2+w0*w0))/(4*pi) |
| | 1384 | |
| | 1385 | s0 <- 4 |
| | 1386 | #s0unit <- s0*(w0+sqrt(2+w0*w0))/(4*pi) |
| | 1387 | s0unit=s0/dt*w0/(2*pi) #(CORRECT) |
| | 1388 | s0log <- as.integer((log2(s0unit)-1)*nvoice+1.5) |
| | 1389 | |
| | 1390 | totnoct <- noctave+as.integer(s0log/nvoice)+1 |
| | 1391 | |
| | 1392 | x <- morlet(N,N/2,sunit,w0) |
| | 1393 | |
| | 1394 | totts.cwt <- Mod(cwt(x,totnoct,nvoice,w0,plot=0)) |
| | 1395 | wt=totts.cwt[,s0log:(s0log+noctave*nvoice)] |
| | 1396 | |
| | 1397 | wt <- wt/max(wt) |
| | 1398 | |
| | 1399 | if (plot==TRUE) plotwt(wt,0,t,s0,noctave,w0,units="",plottitle="Reproducing Kernel") |
| | 1400 | |
| | 1401 | wt |
| | 1402 | |
| | 1403 | } |
| | 1404 | |
| | 1405 | ################### |
| | 1406 | |
| | 1407 | addvalues <- function(nvoice,dnoctave,x,value){ |
| | 1408 | |
| | 1409 | nr <- dim(x)[1] |
| | 1410 | nc <- dim(x)[2] |
| | 1411 | dnc <- nvoice*dnoctave |
| | 1412 | |
| | 1413 | y <- matrix(rep(value,nr*(nc+dnc)),nrow=nr) |
| | 1414 | |
| | 1415 | y[,(dnc+1):(dnc+nc)] <- x |
| | 1416 | |
| | 1417 | y |
| | 1418 | |
| | 1419 | } |
| | 1420 | |
| | 1421 | #################### |
| | 1422 | |
| | 1423 | scalematrix <- function(wt){ |
| | 1424 | |
| | 1425 | # scales matrix, such that the maximal value is one |
| | 1426 | # wt: matrix to be scaled |
| | 1427 | |
| | 1428 | mval <- max(wt,na.rm=TRUE) |
| | 1429 | |
| | 1430 | wt <- wt/mval |
| | 1431 | |
| | 1432 | wt |
| | 1433 | |
| | 1434 | } |
| | 1435 | |
| | 1436 | |
| | 1437 | foldKernel <- function(F,swabs,tw,s,dt){ |
| | 1438 | |
| | 1439 | # folds a matrix (e.g. kernel with smoothing window |
| | 1440 | # F: matrix input |
| | 1441 | # swabs: smooting window width |
| | 1442 | |
| | 1443 | smsF <- smooth.matrix(F,swabs) |
| | 1444 | smsF[is.na(smsF)] <- 0 |
| | 1445 | |
| | 1446 | smtsF <- smooth.time(smsF,tw,dt,s) |
| | 1447 | |
| | 1448 | smtsF[is.na(smtsF)] <- 0 |
| | 1449 | |
| | 1450 | scF <- scalematrix(smtsF) |
| | 1451 | |
| | 1452 | scF |
| | 1453 | |
| | 1454 | } |
| | 1455 | |
| | 1456 | kernelBitmap <- function(c,s0=1,w0=6,noctave=6,nvoice=20,swabs=0,tw=0,dt=1){ |
| | 1457 | # produces kernel bitmap |
| | 1458 | # c: height of contour, that defines area |
| | 1459 | # s0: lowest scale |
| | 1460 | # noctave: number of octaves |
| | 1461 | # nvoice: number of voices per octave |
| | 1462 | # swabs: smoothing window length in scale direction |
| | 1463 | # dt: sampling time |
| | 1464 | |
| | 1465 | s <- s0*2^noctave |
| | 1466 | is <- noctave*nvoice |
| | 1467 | |
| | 1468 | x <- s0*2^(((1:(nvoice*(noctave+2)))-1)/nvoice) |
| | 1469 | |
| | 1470 | t <- max(2000,max(x)*tw*2/dt*1.1) |
| | 1471 | |
| | 1472 | y <- ((0:(2*t))-t)/2*dt |
| | 1473 | |
| | 1474 | X <- matrix(x,ncol=nvoice*(noctave+2),nrow=2*t+1,byrow=T) |
| | 1475 | Y <- matrix(y,ncol=nvoice*(noctave+2),nrow=2*t+1) |
| | 1476 | |
| | 1477 | F <- sqrt(2*s*X/(s*s+X*X))*exp(-0.5*(Y*Y+w0*w0*(X-s)*(X-s))/(s*s+X*X)) |
| | 1478 | |
| | 1479 | F <- foldKernel(F,swabs,tw,x,dt) |
| | 1480 | |
| | 1481 | F[F<c] <- 0 |
| | 1482 | F[F>=c] <- 1 |
| | 1483 | |
| | 1484 | is1 <- 1 |
| | 1485 | is2 <- nvoice*(noctave+1) |
| | 1486 | it1 <- 1 |
| | 1487 | it2 <- 2*t+1 |
| | 1488 | |
| | 1489 | L <- F[1:(2*t+1),is1] |
| | 1490 | |
| | 1491 | while (length(L[L!=0])==0) { |
| | 1492 | is1 <- is1+1 |
| | 1493 | L <- F[1:(2*t+1),is1] |
| | 1494 | } |
| | 1495 | |
| | 1496 | L <- F[1:(2*t+1),is2] |
| | 1497 | |
| | 1498 | while (length(L[L!=0])==0) { |
| | 1499 | is2 <- is2-1 |
| | 1500 | L <- F[1:(2*t+1),is2] |
| | 1501 | } |
| | 1502 | |
| | 1503 | |
| | 1504 | L <- F[it1,1:(nvoice*(noctave+2))] |
| | 1505 | |
| | 1506 | while (length(L[L!=0])==0) { |
| | 1507 | it1 <- it1+1 |
| | 1508 | L <- F[it1,1:(nvoice*(noctave+2))] |
| | 1509 | } |
| | 1510 | |
| | 1511 | L <- F[it2,1:(nvoice*(noctave+2))] |
| | 1512 | |
| | 1513 | while (length(L[L!=0])==0) { |
| | 1514 | it2 <- it2-1 |
| | 1515 | L <- F[it2,1:(nvoice*(noctave+2))] |
| | 1516 | } |
| | 1517 | |
| | 1518 | kernel <- list(bitmap=F[(it1-1):(it2+1),(is1-1):(is2+1)],is=is-is1) |
| | 1519 | |
| | 1520 | kernel |
| | 1521 | |
| | 1522 | } |
| | 1523 | |
| | 1524 | kernelRoot <- function(s0=1,w0=6,a=0,noctave=6,nvoice=20,swabs=0,tw=0,dt=1){ |
| | 1525 | |
| | 1526 | tol <- 0.005 |
| | 1527 | cmin <- 0 |
| | 1528 | cmax <- 1 |
| | 1529 | cntr <- 0.5 |
| | 1530 | da <- a |
| | 1531 | |
| | 1532 | while (abs(da/a)>tol){ |
| | 1533 | |
| | 1534 | da <- kernelArea(cntr,s0,w0,a,noctave,nvoice,swabs,tw,dt) |
| | 1535 | if (da>0){ |
| | 1536 | cmin <- cntr |
| | 1537 | cntr <- mean(c(cntr,cmax)) |
| | 1538 | } |
| | 1539 | if (da<0){ |
| | 1540 | cmax <- cntr |
| | 1541 | cntr <- mean(c(cntr,cmin)) |
| | 1542 | } |
| | 1543 | } |
| | 1544 | |
| | 1545 | cntr |
| | 1546 | |
| | 1547 | } |
| | 1548 | |
| | 1549 | kernelArea <- function(cntr,s0=1,w0=6,a=0,noctave=6,nvoice=20,swabs=0,tw=0,dt=1){ |
| | 1550 | |
| | 1551 | # calulates area of reproducing kernel for smoothed data at scale s0*2^noctave |
| | 1552 | # cntr: height of contour line to define kernel area. This |
| | 1553 | # parameter is to be estimated! |
| | 1554 | # s0: lowest scale |
| | 1555 | # w0: parameter of Morlet Wavelet |
| | 1556 | # a: area offset (only needed, when finding root. Is set to |
| | 1557 | # desired area |
| | 1558 | # noctave: number of octaves |
| | 1559 | # nvoice: number of voices per octave |
| | 1560 | # swabs: smoothing window width in scale direction |
| | 1561 | # dt: sampling time |
| | 1562 | |
| | 1563 | s <- s0*2^noctave |
| | 1564 | |
| | 1565 | x <- s0*2^(((1:(nvoice*(noctave+2)))-1)/nvoice) |
| | 1566 | |
| | 1567 | t <- max(2000,max(x)*tw*2/dt*1.1) |
| | 1568 | |
| | 1569 | y <- ((0:(2*t))-t)/2*dt |
| | 1570 | |
| | 1571 | X <- matrix(x,ncol=nvoice*(noctave+2),nrow=2*t+1,byrow=T) |
| | 1572 | Y <- matrix(y,ncol=nvoice*(noctave+2),nrow=2*t+1) |
| | 1573 | |
| | 1574 | F <- sqrt(2*s*X/(s*s+X*X))*exp(-0.5*(Y*Y+w0*w0*(X-s)*(X-s))/(s*s+X*X)) |
| | 1575 | |
| | 1576 | F <- foldKernel(F,swabs,tw,x,dt) |
| | 1577 | |
| | 1578 | F[F>=cntr] <- 1 |
| | 1579 | F[F<cntr] <- 0 |
| | 1580 | |
| | 1581 | area <- length(F[F==1])-a |
| | 1582 | |
| | 1583 | area |
| | 1584 | |
| | 1585 | } |
| | 1586 | |
| | 1587 | |
| | 1588 | tobin <- function(x){ |
| | 1589 | |
| | 1590 | # sets nonzero values to one |
| | 1591 | |
| | 1592 | y <- x/x |
| | 1593 | y[is.na(y)] <- 0 |
| | 1594 | |
| | 1595 | y |
| | 1596 | |
| | 1597 | } |
| | 1598 | |
| | 1599 | scaleKernel <- function(kernel,l){ |
| | 1600 | |
| | 1601 | # scales kernel length in time direction proportional to scale |
| | 1602 | # kernel: data bitmap of width n |
| | 1603 | # l: new width of kernel |
| | 1604 | |
| | 1605 | n <- nrow(kernel) |
| | 1606 | m <- ncol(kernel) |
| | 1607 | |
| | 1608 | newKernel <- matrix(rep(0,m*l),nrow=l) |
| | 1609 | |
| | 1610 | d <- as.double(n)/as.double(l) |
| | 1611 | |
| | 1612 | for (i in 1:l){ |
| | 1613 | j <- as.integer((i-0.5)*d) |
| | 1614 | if (j==0) j <- 1 |
| | 1615 | newKernel[i,1:m] <- kernel[j,1:m] |
| | 1616 | } |
| | 1617 | |
| | 1618 | newKernel |
| | 1619 | |
| | 1620 | } |
| | 1621 | |
| | 1622 | slope <- function(w0,swabs,tw,nvoice,siglevel,arealsiglevel,type){ |
| | 1623 | |
| | 1624 | sw <- swabs/nvoice |
| | 1625 | |
| | 1626 | if (type==0){ # wavelet spectrum |
| | 1627 | |
| | 1628 | if (tw == 0 & sw == 0 & w0 == 1 *pi) slp <- 5.82518 # w = 18.35831 |
| | 1629 | if (tw == 1.5 & sw == 0 & w0 == 1 *pi) slp <- 24.69852 # w = 14.30709 |
| | 1630 | if (tw == 3 & sw == 0 & w0 == 1 *pi) slp <- 35.48368 # w = 14.72354 |
| | 1631 | if (tw == 0 & sw == 5 & w0 == 1 *pi) slp <- 7.347707 # w = 17.96942 |
| | 1632 | if (tw == 1.5 & sw == 5 & w0 == 1 *pi) slp <- 28.24291 # w = 12.65993 |
| | 1633 | if (tw == 3 & sw == 5 & w0 == 1 *pi) slp <- 51.13723 # w = 10.96359 |
| | 1634 | if (tw == 0 & sw == 10 & w0 == 1 *pi) slp <- 10.47856 # w = 15.5941 |
| | 1635 | if (tw == 1.5 & sw == 10 & w0 == 1 *pi) slp <- 45.07387 # w = 15.29793 |
| | 1636 | if (tw == 3 & sw == 10 & w0 == 1 *pi) slp <- 52.82886 # w = 12.72361 |
| | 1637 | |
| | 1638 | if (tw == 0 & sw == 0 & w0 == 2 *pi) slp <- 8.718912 # w = 17.75933 |
| | 1639 | if (tw == 1.5 & sw == 0 & w0 == 2 *pi) slp <- 11.88006 # w = 15.39648 |
| | 1640 | if (tw == 3 & sw == 0 & w0 == 2 *pi) slp <- 26.55977 # w = 1.064384 |
| | 1641 | if (tw == 0 & sw == 5 & w0 == 2 *pi) slp <- 14.64761 # w = 16.27518 |
| | 1642 | if (tw == 1.5 & sw == 5 & w0 == 2 *pi) slp <- 28.27798 # w = 14.57059 |
| | 1643 | if (tw == 3 & sw == 5 & w0 == 2 *pi) slp <- 63.54121 # w = 12.83778 |
| | 1644 | if (tw == 0 & sw == 10 & w0 == 2 *pi) slp <- 27.78735 # w = 11.95813 |
| | 1645 | if (tw == 1.5 & sw == 10 & w0 == 2 *pi) slp <- 41.27260 # w = 12.03379 |
| | 1646 | if (tw == 3 & sw == 10 & w0 == 2 *pi) slp <- 67.37015 # w = 10.63935 |
| | 1647 | |
| | 1648 | } |
| | 1649 | |
| | 1650 | if (type==1){ # wavelet coherence |
| | 1651 | |
| | 1652 | if (tw==0 & sw==0 & w0==pi) slp <- 999 #not valid |
| | 1653 | if (tw==1.5 & sw==0 & w0==pi) slp <- 1 |
| | 1654 | if (tw==3 & sw==0 & w0==pi) slp <- 1 |
| | 1655 | if (tw==0 & sw==0.5 & w0==pi) slp <- 1 |
| | 1656 | if (tw==1.5 & sw==0.5 & w0==pi) slp <- 1 |
| | 1657 | if (tw==3 & sw==0.5 & w0==pi) slp <- 1 |
| | 1658 | if (tw==0 & sw==1 & w0==pi) slp <- 1 |
| | 1659 | if (tw==1.5 & sw==1 & w0==pi) slp <- 1 |
| | 1660 | if (tw==3 & sw==1 & w0==pi) slp <- 1 |
| | 1661 | |
| | 1662 | if (tw==0 & sw==0 & w0==2*pi) slp <- 999 #not valid |
| | 1663 | if (tw==1.5 & sw==0 & w0==2*pi) slp <- 1 |
| | 1664 | if (tw==3 & sw==0 & w0==2*pi) slp <- 1 |
| | 1665 | if (tw==0 & sw==0.5 & w0==2*pi) slp <- 1 |
| | 1666 | if (tw==1.5 & sw==0.5 & w0==2*pi) slp <- 8.3 |
| | 1667 | if (tw==3 & sw==0.5 & w0==2*pi) slp <- 1 |
| | 1668 | if (tw==0 & sw==1 & w0==2*pi) slp <- 1 |
| | 1669 | if (tw==1.5 & sw==1 & w0==2*pi) slp <- 1 |
| | 1670 | if (tw==3 & sw==1 & w0==2*pi) slp <- 1 |
| | 1671 | |
| | 1672 | if (tw==0 & sw==0 & w0==3*pi) slp <- 999 #not valid |
| | 1673 | if (tw==1.5 & sw==0 & w0==3*pi) slp <- 1 |
| | 1674 | if (tw==3 & sw==0 & w0==3*pi) slp <- 1 |
| | 1675 | if (tw==0 & sw==0.5 & w0==3*pi) slp <- 1 |
| | 1676 | if (tw==1.5 & sw==0.5 & w0==3*pi) slp <- 1 |
| | 1677 | if (tw==3 & sw==0.5 & w0==3*pi) slp <- 1 |
| | 1678 | if (tw==0 & sw==1 & w0==3*pi) slp <- 1 |
| | 1679 | if (tw==1.5 & sw==1 & w0==3*pi) slp <- 1 |
| | 1680 | if (tw==3 & sw==1 & w0==3*pi) slp <- 1 |
| | 1681 | |
| | 1682 | if (tw==0 & sw==0 & w0==4*pi) slp <- 999 #not valid |
| | 1683 | if (tw==1.5 & sw==0 & w0==4*pi) slp <- 1 |
| | 1684 | if (tw==3 & sw==0 & w0==4*pi) slp <- 1 |
| | 1685 | if (tw==0 & sw==0.5 & w0==4*pi) slp <- 1 |
| | 1686 | if (tw==1.5 & sw==0.5 & w0==4*pi) slp <- 1 |
| | 1687 | if (tw==3 & sw==0.5 & w0==4*pi) slp <- 1 |
| | 1688 | if (tw==0 & sw==1 & w0==4*pi) slp <- 1 |
| | 1689 | if (tw==1.5 & sw==1 & w0==4*pi) slp <- 1 |
| | 1690 | if (tw==3 & sw==1 & w0==4*pi) slp <- 1 |
| | 1691 | |
| | 1692 | } |
| | 1693 | |
| | 1694 | cat(paste("# slope ",slp,"\n",sep="")) |
| | 1695 | |
| | 1696 | slp |
| | 1697 | |
| | 1698 | } |
| | 1699 | |
| | 1700 | #============================ |
| | 1701 | #============================ |
| | 1702 | #============================ |
| | 1703 | |
| | 1704 | ## 1. Read auxiliar data files #### |
| | 1705 | |
| | 1706 | identifiants <- read.table("identifs.txt")[ ,1] |
| | 1707 | dates0 <- read.table("datesall.txt")[, 1] |
| | 1708 | dates <- as.character(dates0[grep("2009", dates0)]) |
| | 1709 | rm(dates0) |
| | 1710 | |
| | 1711 | n <- length(identifiants) |
| | 1712 | p <- delta <- length(dates) |
| | 1713 | |
| | 1714 | synchros09 <- t(as.matrix(read.table("~/tmp/2009_synchros200RC.txt"))) |
| | 1715 | #synchros09 <- t(as.matrix(read.table("~/tmp/2009_synchros200-random.txt"))) |
| | 1716 | |
| | 1717 | nas <- which(is.na(synchros09)[, 1]) # some 1/1/2009 are missing |
| | 1718 | synchros09[nas, 1] <- rowMeans(synchros09[nas, 2:4]) #valeurs après 1er janvier |
| | 1719 | |
| | 1720 | #moyenne pondérée pour compléter deux demi-heures manquantes |
| | 1721 | imput09 <- synchros09[, 4180:4181] %*% matrix(c(2/3, 1/3, 1/3, 2/3), 2) |
| | 1722 | synchros09 <- cbind(synchros09[, 1:4180], imput09, synchros09[, 4181:17518]) |
| | 1723 | |
| | 1724 | conso <- synchros09[-201, ]; # series must be on rows |
| | 1725 | n <- nrow(conso) |
| | 1726 | delta <- ncol(conso) |
| | 1727 | |
| | 1728 | rm(synchros09, nas) |
| | 1729 | |
| | 1730 | ## 2. Compute WER distance matrix #### |
| | 1731 | |
| | 1732 | ## _.a CWT -- Filtering the lowest freqs (>6m) #### |
| | 1733 | # nvoice <- 4 |
| | 1734 | # # noctave4 = 2^13 = 8192 half hours ~ 180 days |
| | 1735 | # noctave4 <- adjust.noctave(N = delta, dt = 1, s0 = 2, |
| | 1736 | # tw = 0, noctave = 13) |
| | 1737 | # # 4 here represent 2^5 = 32 half-hours ~ 1 day |
| | 1738 | # scalevector4 <- 2^(4:(noctave4 * nvoice) / nvoice) * 2 |
| | 1739 | # lscvect4 <- length(scalevector4) |
| | 1740 | # lscvect <- lscvect4 # i should clean my code: werFam demands a lscvect |
| | 1741 | |
| | 1742 | |
| | 1743 | #17000 colonnes coeff 1, puis 17000 coeff 2... [non : dans chaque tranche du cube] |
| | 1744 | |
| | 1745 | #TODO: une fonction qui fait lignes 59 à 91 |
| | 1746 | |
| | 1747 | #cube: |
| | 1748 | # Xcwt4 <- toCWT(conso, noctave = noctave4, dt = 1, |
| | 1749 | # scalevector = scalevector4, |
| | 1750 | # lt = delta, smooth = FALSE, |
| | 1751 | # nvoice = nvoice) # observations node with CWT |
| | 1752 | # |
| | 1753 | # #matrix: |
| | 1754 | # ############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect) |
| | 1755 | # #Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1])))) |
| | 1756 | # |
| | 1757 | # #NOTE: delta et lscvect pourraient etre gardés à part (communs) |
| | 1758 | # for(i in 1:n) |
| | 1759 | # Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) ) |
| | 1760 | # |
| | 1761 | # #rm(conso, Xcwt4); gc() |
| | 1762 | # |
| | 1763 | # ## _.b WER^2 distances ######## |
| | 1764 | # Xwer_dist <- matrix(0.0, n, n) |
| | 1765 | # for(i in 1:(n - 1)){ |
| | 1766 | # mat1 <- vect2mat(Xcwt2[i,]) |
| | 1767 | # for(j in (i + 1):n){ |
| | 1768 | # mat2 <- vect2mat(Xcwt2[j,]) |
| | 1769 | # num <- Mod(mat1 * Conj(mat2)) |
| | 1770 | # WX <- Mod(mat1 * Conj(mat1)) |
| | 1771 | # WY <- Mod(mat2 * Conj(mat2)) |
| | 1772 | # smsmnum <- smCWT(num, scalevector = scalevector4) |
| | 1773 | # smsmWX <- smCWT(WX, scalevector = scalevector4) |
| | 1774 | # smsmWY <- smCWT(WY, scalevector = scalevector4) |
| | 1775 | # wer2 <- sum(colSums(smsmnum)^2) / |
| | 1776 | # sum( sum(colSums(smsmWX) * colSums(smsmWY)) ) |
| | 1777 | # Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2)) |
| | 1778 | # Xwer_dist[j, i] <- Xwer_dist[i, j] |
| | 1779 | # } |
| | 1780 | # } |
| | 1781 | # diag(Xwer_dist) <- numeric(n) |
| | 1782 | # |
| | 1783 | # save(Xwer_dist, file = "../res/2009_synchros200WER.Rdata") |
| | 1784 | # save(Xwer_dist, file = "../res/2009_synchros200-randomWER.Rdata") |
| | 1785 | |
| | 1786 | load("../res/2009_synchros200WER.Rdata") |
| | 1787 | #load("../res/2009_synchros200-randomWER.Rdata") |
| | 1788 | |
| | 1789 | ## 3. Cluster using WER distance matrix #### |
| | 1790 | |
| | 1791 | #hc <- hclust(as.dist(Xwer_dist), method = "ward.D") |
| | 1792 | #plot(hc) |
| | 1793 | # |
| | 1794 | # #clust <- cutree(hc, 2) |
| | 1795 | # |
| | 1796 | for(K in 2:30){ |
| | 1797 | #K <- 3 |
| | 1798 | #pamfit <- pam(tdata[-201, ci$selectv], k = K) |
| | 1799 | pamfit <- pam(as.dist(Xwer_dist), k = K, diss = TRUE) |
| | 1800 | |
| | 1801 | #table(pamfit$clustering) |
| | 1802 | |
| | 1803 | SC <- matrix(0, ncol = p, nrow = K) |
| | 1804 | |
| | 1805 | clustfactor <- pamfit$clustering |
| | 1806 | # for(k in 1:K){ |
| | 1807 | # clustk <- which(clustfactor == k) |
| | 1808 | # if(length(clustk) > 0) { |
| | 1809 | # if(length(clustk) > 1) { |
| | 1810 | # SCk <- colSums(synchros09[which(clustfactor == k), ]) |
| | 1811 | # } else { |
| | 1812 | # SCk <- synchros09[which(clustfactor == k), ] |
| | 1813 | # } |
| | 1814 | # SC[k, ] <- SC[k, ] + SCk |
| | 1815 | # rm(SCk) |
| | 1816 | # } |
| | 1817 | #} |
| | 1818 | |
| | 1819 | #write.table(clustfactor, file = paste0("~/tmp/clustfactorRC", K, ".txt")) |
| | 1820 | #write.table(clustfactor, file = "~/tmp/clustfactor3.txt") |
| | 1821 | #write.table(clustfactor, file = paste0("~/tmp/clustfactorWER", K, ".txt")) |
| | 1822 | write.table(clustfactor, file = paste0("~/tmp/clustfactor-randomWER", K, ".txt")) |
| | 1823 | } |
| | 1824 | # |
| | 1825 | # # Plots |
| | 1826 | # layout(1) |
| | 1827 | # matplot(t(SC)[48*10 + 1:(48*30), ], type = 'l', ylab = '',col = 1:3, lty = 1) |
| | 1828 | # matplot(t(SC)[48*100 + 1:(48*30), ], type = 'l', ylab = '', col = 1:3, lty = 1) |
| | 1829 | # |
| | 1830 | # |
| | 1831 | # |
projects / epclust.git / blame_incremental