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