X-Git-Url: https://git.auder.net/assets/current/gitweb.css?a=blobdiff_plain;f=old_C_code%2Fstage2_UNFINISHED%2Fsrc%2Funused%2Fsowas-superseded.r;fp=old_C_code%2Fstage2_UNFINISHED%2Fsrc%2Funused%2Fsowas-superseded.r;h=0000000000000000000000000000000000000000;hb=62deb4244895a20a35397dfb062f0b9fe94c5012;hp=dca33536ce316493b67ff989bd2f7064e4d45580;hpb=3eef8d3df59ded9a281cff51f79fe824198a7427;p=epclust.git diff --git a/old_C_code/stage2_UNFINISHED/src/unused/sowas-superseded.r b/old_C_code/stage2_UNFINISHED/src/unused/sowas-superseded.r deleted file mode 100644 index dca3353..0000000 --- a/old_C_code/stage2_UNFINISHED/src/unused/sowas-superseded.r +++ /dev/null @@ -1,1563 +0,0 @@ -########################################################### -# CONTINUOUS WAVELET TRANSFORMATION OF A TIME SERIES OBJECT -########################################################### - -cwt.ts <- function(ts,s0,noctave=5,nvoice=10,w0=2*pi){ - - if (class(ts)!="ts"){ - - cat("# This function needs a time series object as input. You may construct this by using the function ts(data,start,deltat). Try '?ts' for help.\n") - - } - else{ - - t=time(ts) - dt=t[2]-t[1] - - s0unit=s0/dt*w0/(2*pi) - s0log=as.integer((log2(s0unit)-1)*nvoice+1.5) - - if (s0log<1){ - cat(paste("# s0unit = ",s0unit,"\n",sep="")) - cat(paste("# s0log = ",s0log,"\n",sep="")) - cat("# s0 too small for w0! \n") - } - totnoct=noctave+as.integer(s0log/nvoice)+1 - - totts.cwt=cwt(ts,totnoct,nvoice,w0,plot=0) - - ts.cwt=totts.cwt[,s0log:(s0log+noctave*nvoice)] - - #Normalization - sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0) - smat <- matrix(rep(sqs,length(t)),nrow=length(t),byrow=TRUE) - - ts.cwt*smat - - } - -} - -##################################### -# WSP -##################################### - -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){ - - if (class(ts)!="ts"){ - - cat("# This function needs a time series object as input. You may construct this by using the function ts(data,start,deltat). Try '?ts' for help.\n") - - } - else{ - - if (sw!=0 & swabs==0) - swabs <- as.integer(sw*nvoice) - if (swabs!=0 & sw==0) - sw <- swabs/nvoice - - sllist <- checkarealsiglevel(sw,tw,w0,arealsiglevel,siglevel,0) - arealsiglevel <- sllist$arealsiglevel - siglevel <- sllist$siglevel - - at <- NULL - - t <- time(ts) - dt <- deltat(ts) - s0rem <- s0 - s0 <- adjust.s0(s0,dt) - dnoctave <- as.integer(log(s0/s0rem-0.000000000001)/log(2))+1 - - noctave <- adjust.noctave(length(ts),dt,s0,tw,noctave) - scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0 - tsanom <- ts-mean(ts) - - #WAVELET TRANSFORMATION - ts.cwt <- cwt.ts(tsanom,s0,noctave,nvoice,w0) - - #SMOOTHING - wsp <- smooth.matrix(Mod(ts.cwt)^2,swabs) - smwsp <- smooth.time(wsp,tw,dt,scalevector) - - #POINTWISE SIGNIFICANCE TEST - if (is.null(critval)==FALSE){ # is critval empty? - if (dim(critval)[2]!=dim(smwsp)[2]){ # is critval of the wrong dimension? - if (siglevel[1]!=0 & nreal!=0) critval <- - criticalvaluesWSP(tsanom,s0,noctave,nvoice,w0,swabs,tw,siglevel,nreal) - #critval is of wrong dimension and siglevel and nreal are given - else { - critval <- NULL # no test possible - arealsiglevel <- 0 - cat("# dimension of critval is wrong \n") - cat("# areawise test only possible with pointwise test \n") - } - } - } - else{ # critval is empty, but nreal or siglevel is given - if (siglevel[1]!=0 & nreal!=0) critval <- - criticalvaluesWSP(tsanom,s0,noctave,nvoice,w0,swabs,tw,siglevel,nreal) - else { - critval <- NULL - arealsiglevel <- 0 - cat("# areawise test only possible with pointwise test \n") - } - } - - #AREAL SIGNIFICANCE TEST - if (arealsiglevel!=0){ - v <- critval[1,] - v[v==0] <- 9999 - cvmat <- matrix(rep(v,length(t)),nrow=length(t),byrow=TRUE) - atest <- arealtest(smwsp/cvmat,dt,s0,noctave,nvoice,w0,swabs,tw,siglevel,arealsiglevel,kernel,0) - at <- atest$at - kernel <- atest$kernel - } - - if (s0rem1) cv[i] <- 1 - - cat(paste("# significance testing, cv=",cv[i]," \n",sep="")) - - } - - cv - -} - -############################# -# AREAWISE SIGNIFICANCE TEST -############################# - -slide <- function(data,kernellist,s0,noctave,nvoice,cv){ - - # slides kernel over data set - #---------------------------- - # data: matrix containing data - # kernellist: matrix containing kernel - # s0: lowest scale - # noctave: number of octaves - # nvoice: number of voices per octave - # cv: critical value, all values higher are set to one - - #Time: (rows) n,i the kernel is to be scaled in this direction - #Scale: (cols) m,j - - data[data1){ - linewidth <- 1 - if (sigplot==3) - linewidth <- 5 - - contour(x=t,z=at,levels=0.5,add=TRUE,drawlabels=FALSE,axes=FALSE,lwd=linewidth,col="black") - } - -} - - -plotcv <- function(t,wt,cv){ - - if (length(dim(cv))==0) - - contour(x=t,z=wt,levels=c(cv),drawlabels=FALSE,axes=FALSE,add=TRUE,col="black",lwd=1) - - else{ - - for(i in 1:nrow(cv)){ - - v <- cv[i,] - v[v==0] <- 9999 - m <- matrix(rep(v,length(t)),nrow=length(t),byrow=TRUE) - - contour(x=t,z=wt/m,levels=1,drawlabels=FALSE,axes=FALSE,add=TRUE,col="black",lwd=1) - - } - - } - -} - -plotcoi <- function(t,s0,noctave,w0){ - - tv <- as.vector(t) - tvl <- tv[tv-tv[1]<(tv[length(tv)]-tv[1])/2] - tvr <- tv[tv-tv[1]>=(tv[length(tv)]-tv[1])/2] - - lines(tvl,log2(((tvl-tv[1])*4*pi/((w0+sqrt(2+w0*w0))*sqrt(2)))/s0)/noctave,col="black") - lines(tvr,log2(((tv[length(tv)]-tvr)*4*pi/((w0+sqrt(2+w0*w0))*sqrt(2)))/s0)/noctave,col="black") - -} - -plotmarks <- function(t,s0,noctave,markt,marks){ - - if (markt[1]!=-999){ - - for (i in 1:length(markt)){ - lines(c(markt[i],markt[i]),c(0,1),lty="dotted") - } - - } - - if (marks[1]!=-999){ - - for (i in 1:length(marks)){ - lines(c(t[1],t[length(t)]),c(log2(marks[i]/s0)/noctave,log2(marks[i]/s0)/noctave),lty="dotted") - } - - } - -} - - -##################### -# PLOT.WT -##################### - -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){ - - 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) - -} - -plotwt <- - 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){ - - if (is.null(phs)) phase <- FALSE - - mgpv <- c(3,1,0) - if (labsc>1) mgpv[1] <- 3-(labsc-1.5) - - ncolors <- 64 - colors <- wtcolors(ncolors) - if (color==FALSE) colors <- gray((0:ncolors)/ncolors*0.7+0.3) - - rangev <- (0:(ncolors-1))/(ncolors-1) - rangebar <- matrix(rangev,nrow=2,ncol=64,byrow=TRUE) - - s <- 2^(0:(noctave))*s0 - sn <- (0:(noctave))/noctave - - deltat <- (t[length(t)]-t[1])/(length(t)-1) - cut <- FALSE - if ((!is.null(t1)) | (!is.null(t2))){ - if (t1=t[1] & t2<=t[length(t)]){ - - cut <- TRUE - - i1 <- (t1-t[1])/deltat+1 - i2 <- (t2-t[1])/deltat+1 - - t <- t[t>=t1 & t<=t2] - - wt <- wt[i1:i2,] - if (phase) phs <- phs[i1:i2,] - if (length(at)>1) at <- at[i1:i2,] - - } - } - - #---------------- - # Setting layout - #---------------- - - if (device=="ps" && split==FALSE){ - - file <- paste(file,".eps",sep="") - - postscript(file,onefile=TRUE,horizontal=TRUE,paper="special",width=pwidth,height=pheight) - cat(paste("# writing",file," \n")) - - } - - if (device=="ps" && split==TRUE){ - - file1 <- paste(file,".mod.eps",sep="") - - postscript(file1,onefile=TRUE,horizontal=TRUE,paper="special",width=pwidth,height=pheight) - cat(paste("# writing",file1," \n")) - - } - - if (phase==TRUE && split==FALSE) nlo <- layout(matrix(c(1,2,3,4),2,2,byrow=TRUE),widths=c(4,1)) - else nlo <- layout(matrix(c(1,2),ncol=2,byrow=TRUE),widths=c(4,1)) - - - #----------------- - # Plotting modulus - #----------------- - - if (logscale){ - if (units=="") - image(x=t,z=log10(wt),col = colors,axes=FALSE,xlab="Time",ylab="Scale",frame.plot=TRUE,cex.lab=labsc,mgp=mgpv) - else - 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) - } - else{ - if (units=="") - image(x=t,z=wt,col = colors,axes=FALSE,xlab="Time",ylab="Scale",frame.plot=TRUE,cex.lab=labsc,mgp=mgpv) - else - 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) - } - - text(t[1+as.integer(length(t)*0.1)],0.8,labtext,cex=2) - - if (sigplot==1 | sigplot==3){ - if (is.null(cv)==FALSE){ #Critical values - if (cv[1]!=0 & cv[1]!=-1) plotcv(t,wt,cv) - } - } - - if (sigplot>1) plotat(t,wt,at,sigplot) - - if (!cut) plotcoi(t,s0,noctave,w0) #cone of influence - plotmarks(t,s0,noctave,markt,marks) #additional guiding lines - - if (is.null(xax)) - axis(side=1,at=axTicks(1),cex.axis=labsc) - else - if (is.null(xlab)) - axis(side=1,xax,labels=as.character(xax),cex.axis=labsc) - else - axis(side=1,xax,labels=xlab,cex.axis=labsc) - - if (is.null(yax)) - axis(side=2,sn,labels=as.character(s),cex.axis=labsc) - else - if (is.null(ylab)) - axis(side=2,yax,labels=as.character(yax),cex.axis=labsc) - else - axis(side=2,yax,labels=ylab,cex.axis=labsc) - - title(main=plottitle,cex=labsc) - - image(z=rangebar,axes=FALSE,col=colors,frame.plot=TRUE,cex.lab=labsc,mgp=mgpv) - - if (is.null(cv)==FALSE){ - if (length(dim(cv))==0){ - for (i in 1:length(cv)) - if (cv[i]!=0) lines(c(-1,2),c(cv[i],cv[i])) - } - } - - if (!logscale) - axis(side=2,(0:5)/5,labels=c("0","","","","","1"),cex.axis=labsc) - else{ - labelv <- substr(as.character(c(0:5)*(max(log10(wt),na.rm=TRUE)-min(log10(wt),na.rm=TRUE))/5),1,4) - axis(side=2,(0:5)/5,labels=labelv,cex.axis=labsc) - } - - - #----------------- - # Plotting phase - #----------------- - if (phase==TRUE){ - - if (device=="ps" && split==TRUE){ - - dev.off() - - file2 <- paste(file,".phs.eps",sep="") - - postscript(file2,onefile=TRUE,horizontal=TRUE,paper="special",width=10,height=5) - cat(paste("# writing",file2," \n")) - - } - - colors <- rainbow(ncolors) - if (color==FALSE){ - dummy <- gray((0:ncolors)/ncolors) - colors[1:(ncolors/2)] <- dummy[(ncolors/2+1):ncolors] - colors[(ncolors/2+1):ncolors] <- dummy[1:(ncolors/2)] - } - - if (units=="") - image(x=t,z=phs,col=colors,axes=FALSE,xlab="Time",ylab="Scale",frame.plot=TRUE,cex.lab=labsc,mgp=mgpv) - else - 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) - - if (is.null(cv)==FALSE) plotcv(t,wt,cv) - if (sigplot>1) plotat(t,wt,at,sigplot) - if (!cut) plotcoi(t,s0,noctave,w0) - plotmarks(t,s0,noctave,markt,marks) - - if (is.null(xax)) - axis(side=1,at=axTicks(1),cex.axis=labsc) - else - if (is.null(xlab)) - axis(side=1,xax,labels=as.character(xax),cex.axis=labsc) - else - axis(side=1,xax,labels=xlab,cex.axis=labsc) - - if (is.null(yax)) - axis(side=2,sn,labels=as.character(s),cex.axis=labsc) - else - if (is.null(ylab)) - axis(side=2,yax,labels=as.character(yax),cex.axis=labsc) - else - axis(side=2,yax,labels=ylab,cex.axis=labsc) - - - title(main="Phase") - - image(z=rangebar,axes=FALSE,col=colors,frame.plot=TRUE,cex.lab=labsc,mgp=mgpv) - axis(side=2,(0:4)/4,labels=c("-PI","","0","","PI"),cex.axis=labsc) - - } - - if (device=="ps") dev.off() - - } - -############################## -# Surrogates -############################## - -createwavesurrogates <- function(nsurr=1,wt=1,n,dt=1,s0=1,noctave=5,nvoice=10,w0=2*pi){ - - surrmat <- matrix(rep(0,n*nsurr),ncol=nsurr) - - for (i in 1:nsurr){ - - cat(paste("# Creating realization ",i,"\n",sep="")) - - x <- rnorm(n) - xts <- ts(x,start=0,deltat=dt) - - xwt <- cwt.ts(xts,s0,noctave,nvoice,w0) - wtsurr <- wt*xwt - - surri <- as.vector(invmorlet(wtsurr,0,dt,s0,noctave,nvoice,w0)) - - surrmat[,i] <- Re(surri) - - } - - surrmat - -} - -surrwave <- function(x,...) - UseMethod("surrwave") - -surrwave.wt <- function(wt,nsurr=1,spec=TRUE){ - - n <- length(wt$time) - t0 <- wt$time[1] - dt <- (wt$time[13]-t0)/12 - s0 <- wt$s0 - noctave <- wt$noctave - nvoice <- wt$nvoice - w0 <- wt$w0 - - wt <- wt$modulus - if (spec==TRUE) - wt <- sqrt(wt) - - surrmat <- createwavesurrogates(nsurr,wt,n,dt,s0,noctave,nvoice,w0) - - surrts <- ts(surrmat,start=t0,deltat=dt) - - surrts - -} - -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){ - - if (sw!=0 & swabs==0) - swabs <- as.integer(sw*nvoice) - - scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0 - - if ((noctave*nvoice+1)!=dim(mat)[2]) - cat("# ERROR! nscales unequal noctave*nvoice+1 ! \n") - - n <- dim(mat)[1] - - if (spec==FALSE) - mat <- Mod(mat) - else - mat <- sqrt(Mod(mat)) - - wtsm <- smooth.matrix(mat,swabs) - wt <- smooth.time(wtsm,tw,dt,scalevector) - - surrmat <- createwavesurrogates(nsurr,wt,n,dt,s0,noctave,nvoice,w0) - - surrts <- ts(surrmat,start=t0,deltat=dt) - - surrts - -} - -surrwave.character <- - 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){ - - if (sw!=0 & swabs==0) - swabs <- as.integer(sw*nvoice) - - scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0 - - if (transpose==FALSE) - mat <- matrix(scan(file,comment.char="#"),ncol=nvoice*noctave+1,byrow=TRUE) - else - mat <- matrix(scan(file,comment.char="#"),ncol=nvoice*noctave+1,byrow=FALSE) - - if ((noctave*nvoice+1)!=dim(mat)[2]) - cat("# ERROR! nscales unequal noctave*nvoice+1 ! \n") - - n <- dim(mat)[1] - - if (spec==FALSE) - mat <- Mod(mat) - else - mat <- sqrt(Mod(mat)) - - wtsm <- smooth.matrix(mat,swabs) - wt <- smooth.time(wtsm,tw,dt,scalevector) - - surrmat <- createwavesurrogates(nsurr,wt,n,dt,s0,noctave,nvoice,w0) - - surrts <- ts(surrmat,start=t0,deltat=dt) - - surrts - - } - -surrwave.ts <- function(ts,nsurr=1,s0=1,noctave=5,nvoice=10,w0=2*pi,sw=0,tw=0,swabs=0){ - - n <- length(ts) - t0 <- time(ts)[1] - dt <- deltat(ts) - if (sw!=0 & swabs==0) - swabs <- as.integer(sw*nvoice) - - scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0 - - wt <- Mod(cwt.ts(ts,s0,noctave,nvoice,w0)) - - wtsm <- smooth.matrix(wt,swabs) - wt <- smooth.time(wtsm,tw,dt,scalevector) - - surrmat <- createwavesurrogates(nsurr,wt,n,dt,s0,noctave,nvoice,w0) - - surrts <- ts(surrmat,start=t0,deltat=dt) - - surrts - -} - -invmorlet <- function(wt,t0=0,dt=1,s0=1,noctave=5,nvoice=10,w0=2*pi){ - - if ((noctave*nvoice+1)!=dim(wt)[2]) - cat("# ERROR! nscales unequal noctave*nvoice+1 ! \n") - - n <- dim(wt)[1] - - wt[is.na(wt)] <- 0 - - tsRe <- rep(0,n) - tsIm <- rep(0,n) - - wtRe <- t(Re(wt)) - wtIm <- t(Im(wt)) - - z <- .C("invmorlet", - as.double(as.vector(wtRe)), - as.double(as.vector(wtIm)), - as.integer(n), - as.double(dt), - as.double(s0), - as.integer(noctave), - as.integer(nvoice), - as.double(w0), - tsRetmp = as.double(tsRe), - tsImtmp = as.double(tsIm), - PACKAGE="sowas") - - invvec=complex(real=z$tsRetmp,imaginary=z$tsImtmp) - invts <- ts(data=invvec,start=t0,deltat=dt) - - invts - -} - -################################# -# INPUT / OUTPUT -################################# - -readmatrix <- function(file,M){ - - A <- matrix(scan(file,comment.char="#"),ncol=M,byrow=TRUE) - - A - -} - - -readts <- function(file){ - - A <- matrix(scan(file,comment.char="#"),ncol=2,byrow=TRUE) - - Adum <- A - - Adum[is.na(Adum)] <- 0 - - t <- Adum %*% c(1,0) - x <- A %*% c(0,1) - - N=length(t) - f=1/(t[13]-t[1])*12 - - if ((f>11) && (f<13)) f <- 12 - - timeseries<-ts(data=x,start=t[1],frequency=f) - - timeseries - -} - -writematrix <- function(file,data,header="# R Matrix"){ - - write(header,file) - write(t(data),file,ncol(data),append=TRUE) - -} - -############################ -# HELP FUNCTIONS -############################ - -smooth.matrix <- function(wt,swabs){ - - if (swabs != 0) - smwt <- t(filter(t(wt),rep(1,2*swabs+1)/(2*swabs+1))) - else - smwt <- wt - - smwt - -} - -smooth.time <- function(wt,tw,dt,scalevector){ - - smwt <- wt - - if (tw != 0){ - for (i in 1:length(scalevector)){ - - twi <- as.integer(scalevector[i]*tw/dt) - smwt[,i] <- filter(wt[,i],rep(1,2*twi+1)/(2*twi+1)) - - } - } - - smwt - -} - - -adjust.noctave <- function(N,dt,s0,tw,noctave){ - - if (tw>0){ - dumno <- as.integer((log(N*dt)-log(2*tw*s0))/log(2)) - if (dumno=c] <- 1 - - is1 <- 1 - is2 <- nvoice*(noctave+1) - it1 <- 1 - it2 <- 2*t+1 - - L <- F[1:(2*t+1),is1] - - while (length(L[L!=0])==0) { - is1 <- is1+1 - L <- F[1:(2*t+1),is1] - } - - L <- F[1:(2*t+1),is2] - - while (length(L[L!=0])==0) { - is2 <- is2-1 - L <- F[1:(2*t+1),is2] - } - - - L <- F[it1,1:(nvoice*(noctave+2))] - - while (length(L[L!=0])==0) { - it1 <- it1+1 - L <- F[it1,1:(nvoice*(noctave+2))] - } - - L <- F[it2,1:(nvoice*(noctave+2))] - - while (length(L[L!=0])==0) { - it2 <- it2-1 - L <- F[it2,1:(nvoice*(noctave+2))] - } - - kernel <- list(bitmap=F[(it1-1):(it2+1),(is1-1):(is2+1)],is=is-is1) - - kernel - -} - -kernelRoot <- function(s0=1,w0=6,a=0,noctave=6,nvoice=20,swabs=0,tw=0,dt=1){ - - tol <- 0.005 - cmin <- 0 - cmax <- 1 - cntr <- 0.5 - da <- a - - while (abs(da/a)>tol){ - - da <- kernelArea(cntr,s0,w0,a,noctave,nvoice,swabs,tw,dt) - if (da>0){ - cmin <- cntr - cntr <- mean(c(cntr,cmax)) - } - if (da<0){ - cmax <- cntr - cntr <- mean(c(cntr,cmin)) - } - } - - cntr - -} - -kernelArea <- function(cntr,s0=1,w0=6,a=0,noctave=6,nvoice=20,swabs=0,tw=0,dt=1){ - - # calulates area of reproducing kernel for smoothed data at scale s0*2^noctave - # cntr: height of contour line to define kernel area. This - # parameter is to be estimated! - # s0: lowest scale - # w0: parameter of Morlet Wavelet - # a: area offset (only needed, when finding root. Is set to - # desired area - # noctave: number of octaves - # nvoice: number of voices per octave - # swabs: smoothing window width in scale direction - # dt: sampling time - - s <- s0*2^noctave - - x <- s0*2^(((1:(nvoice*(noctave+2)))-1)/nvoice) - - t <- max(2000,max(x)*tw*2/dt*1.1) - - y <- ((0:(2*t))-t)/2*dt - - X <- matrix(x,ncol=nvoice*(noctave+2),nrow=2*t+1,byrow=T) - Y <- matrix(y,ncol=nvoice*(noctave+2),nrow=2*t+1) - - F <- sqrt(2*s*X/(s*s+X*X))*exp(-0.5*(Y*Y+w0*w0*(X-s)*(X-s))/(s*s+X*X)) - - F <- foldKernel(F,swabs,tw,x,dt) - - F[F>=cntr] <- 1 - F[F