Commit | Line | Data |
---|---|---|
ad642dc6 BA |
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 | } |