projects / epclust.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
complete first draft of package
[epclust.git] / old_C_code / stage2_UNFINISHED / src / unused / sowas-superseded.r
CommitLineData
ad642dc6
BA		 1###########################################################
2# CONTINUOUS WAVELET TRANSFORMATION OF A TIME SERIES OBJECT
3###########################################################
4
5cwt.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
45wsp <- 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
145wco <-
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
264wcs <- 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
343rawWSP <- 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
358criticalvaluesWSP <- 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
423criticalvaluesWCO <- 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
469slide <- 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
528arealtest <- 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
555plotat <- 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
568plotcv <- 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
590plotcoi <- 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
601plotmarks <- 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
626plot.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
632plotwt <-
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
825createwavesurrogates <- 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
849surrwave <- function(x,...)
850 UseMethod("surrwave")
851
852surrwave.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
874surrwave.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
902surrwave.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
936surrwave.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
959invmorlet <- 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
998readmatrix <- function(file,M){
999
1000 A <- matrix(scan(file,comment.char="#"),ncol=M,byrow=TRUE)
1001
1002 A
1003
1004}
1005
1006
1007readts <- 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
1029writematrix <- 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
1040smooth.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
1051smooth.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
1069adjust.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
1083adjust.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
1094adjust.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
1108checkarealsiglevel <- 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
1159as.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
1203wtcolors <- 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
1222createwgn <- function(N,sig,dt){
1223
1224 timeseries<-ts(data=rnorm(N,0,sig),start=0,deltat=dt)
1225
1226 timeseries
1227
1228}
1229
1230
1231createar <- 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
1244rk <- 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
1272addvalues <- 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
1288scalematrix <- 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
1302foldKernel <- 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
1321kernelBitmap <- 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
1389kernelRoot <- 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
1414kernelArea <- 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
1453tobin <- 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
1464scaleKernel <- 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
1487slope <- 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}
Electricity Power Curves Clustering