################################################### ### chunk number 1: loadpacks ################################################### library(wavethresh) library(Rwave) library(ts) library(Milan)

################################################### ### chunk number 2: plotwavelets ################################################### op = par(mfrow=c(5,2), oma = c(0,0, 4, 0), mgp = c(1.2, .8, 0), mar = .1 + c(4,4, .5,1)) for(fn in 1:10) { draw.default(filter.number= fn, col = "blue", main = NULL, xlab= "") abline(h=0,v=0, lty=3, lwd=.5, col = "gray") } mtext(paste("draw.default(*, family = ",formals(draw.default)$family[[2]],")"), side = 3, line = 1, outer = TRUE, cex = par("cex.main"), font = par("font.main")) par(op)

################################################### ### chunk number 3: Example1 ################################################### test.data1=example.1()$y ts.plot(test.data1,ylab="y") title("A function with a jump")

################################################### ### chunk number 4: DWT1 ################################################### wds1 <- wd(test.data1) plot(wds1)

################################################### ### chunk number 5: tree1 ################################################### par(mfrow=c(4,2)) plot.ts(accessC(wds1, level=8)) plot.ts(accessD(wds1, level=8)) plot.ts(accessC(wds1, level=7)) plot.ts(accessD(wds1, level=7)) plot.ts(accessC(wds1, level=6)) plot.ts(accessD(wds1, level=6)) plot.ts(accessC(wds1, level=5)) plot.ts(accessD(wds1, level=5))

################################################### ### chunk number 6: Example2 ################################################### par(mfrow=c(1,1)) test.data2<-achirp(512)$y ts.plot(test.data2,ylab="y") title("A linear chirp")

################################################### ### chunk number 7: DWT2 ################################################### wds2 <- wd(test.data2) plot(wds2)

################################################### ### chunk number 8: tree2 ################################################### par(mfrow=c(4,2)) plot.ts(accessC(wds2, level=8)) plot.ts(accessD(wds2, level=8)) plot.ts(accessC(wds2, level=7)) plot.ts(accessD(wds2, level=7)) plot.ts(accessC(wds2, level=6)) plot.ts(accessD(wds2, level=6)) plot.ts(accessC(wds2, level=5)) plot.ts(accessD(wds2, level=5))

################################################### ### chunk number 9: CWT ################################################### par(mfrow=c(3,1)) plot.ts(test.data2) title("A chirp") cwtchirp<-cwt(test.data2,5,12) image(Arg(cwtchirp))

################################################### ### chunk number 10: Approx1 ################################################### levels<-3:(wds1$nlevels - 1) nthresh <- length(levels) d <- NULL dz <- 0 for (i in 1:nthresh) { d <- accessD(wds1, level = levels[i]) d[abs(d) <= 0.06] = 0 dz=dz+sum(d==0) cat("Level: ", levels[i], " there are ", sum(d == 0), " zeroes\n") wds2 <- putD(wds1, level = levels[i], v = d) } cat("there are in total", dz, " zeroes\n")

################################################### ### chunk number 11: Approx2 ################################################### par(mfrow=c(2,1)) ts.plot(wr(wds2)) # look at the difference ts.plot(test.data1-wr(wds2))

################################################### ### chunk number 12: Denoise1 ################################################### par(mfcol = c(3, 2), mar = c(2, 2, 2, 1), mgp = c(5, 0.4, 0))

x <- (1:128)/128 f1 <- 2 sin(2 pi * x) f2 <- c(rep(0, 32), rep(4, 64), rep(0, 32)) y1 <- f1 + rnorm(128) # performing the dwt of y1 y1wd <- wd(y1, filter.number = 5) plot(x, y1, cex = 0.55) lines(x, wr(threshold(y1wd, policy = "manual", type = "hard", levels = 0:6, value = 1.75))) mtext("Threshold = 1.75", side = 3, line = 0.1, cex = 0.6) plot(x, y1, cex = 0.55) lines(x, wr(threshold(y1wd, policy = "manual", type = "hard", levels = 0:6, value = 2.75))) mtext("Threshold = 2.75", side = 3, line = 0.1, cex = 0.6) plot(x, y1, cex = 0.55) lines(x, wr(threshold(y1wd, policy = "manual", type = "hard", levels = 0:6, value = 3.75))) mtext("Threshold = 3.75", side = 3, line = 0.1, cex = 0.6)

y2 <- f2 + rnorm(128) y2wd <- wd(y2, filter.number = 5) plot(x, y2, cex = 0.55) lines(x, wr(threshold(y2wd, policy = "manual", type = "hard", levels = 0:6, value = 1.75))) mtext("Threshold = 1.75", side = 3, line = 0.1, cex = 0.6) plot(x, y2, cex = 0.55) lines(x, wr(threshold(y2wd, policy = "manual", type = "hard", levels = 0:6, value = 2.75))) mtext("Threshold = 2.75", side = 3, line = 0.1, cex = 0.6) plot(x, y2, cex = 0.55) lines(x, wr(threshold(y2wd, policy = "manual", type = "hard", levels = 0:6, value = 3.75))) mtext("Threshold = 3.75", side = 3, line = 0.1, cex = 0.6)

################################################### ### chunk number 13: Denoise2 ################################################### par(mfrow = c(2, 2), mar = c(1.5, 1.5, 1.5, 0.5), mgp = c(5, 0.4, 0)) x <- (1:2048)/2048 f <- fblocks(x) ssig <- sqrt(var(f)) f <- ((f - mean(f)) * 5)/ssig fnoise <- f + rnorm(2048) fwd <- wd(fnoise, filter.number = 5) plot(x, f, type = "l") mtext(side = 3, line = 0.1, "True function") plot(x, fnoise, type = "l") mtext(side = 3, line = 0.1, "Noisy observations") fwdsure <- threshsure(fwd) plot(x, wr(fwdsure), type = "l") mtext(side = 3, line = 0.1, "SURE estimator") fwdhyb <- threshhyb(fwd, lowlev = 2, seed = 0) plot(x, wr(fwdhyb), type = "l") mtext(side = 3, line = 0.1, "SURE hybrid estimator")

################################################### ### chunk number 14: TI1 ################################################### # k1 and k2 are translation amounts k1 <- 32 k2 <- 42 par(mfrow = c(3, 1), mar = c(1.5, 1.5, 1.5, 0), mgp = c(5, 0.4, 0)) x <- (1:128)/128 unshft <- c(rep(1, 64), rep(-1, 64)) noise <- rnorm(128) # ORIGINAL ESTIMATOR y <- 2 * unshft + noise plot(x, y, type = "l") yest <- wr(threshold(wd(y, filter.number = 1), levels = 0:6)) lines(x, yest) # FIRST TRANSLATED ESTIMATOR if(k1 == 0) y1 <- y else y1 <- c(y[(1 + k1):128], y[1:k1]) if(k1 == 0) shft1 <- unshft else shft2 <- c(unshft[(1 + k1):128], unshft[1:k1]) plot(x, y1, type = "l") y1est <- wr(threshold(wd(y1, filter.number = 1), levels = 0:6)) lines(x, y1est) # SECOND TRANSLATED ESTIMATOR if(k2 == 0) y2 <- y else y2 <- c(y[(1 + k2):128], y[1:k2]) if(k2 == 0) shft2 <- unshft else shft2 <- c(unshft[(1 + k2):128], unshft[1:k2]) plot(x, y2, type = "l") y2est <- wr(threshold(wd(y2, filter.number = 1), levels = 0:6)) lines(x, y2est)

################################################### ### chunk number 15: TI2 ################################################### par(mfrow=c(1,1)) par(mar = c(1.5, 1.5, 1.5, 0.5), mgp = c(5, 0.4, 0)) x <- (1:128)/128 y <- 2 * c(rep(1, 64), rep(-1, 64)) + rnorm(128) plot(x, y, type = "l") yest <- spincyclen(y, filter.number = 1) lines(x, yest) yest1 <- spincyclen(y, filter.number = 3) lines(x, yest1,col=red)