library(waveslim) ### There is a function in waveslim called universal.thresh, ### but it only allows for hard and soft thresholding. The ### following adaptation handles hard, soft and mid thresholding, ### and is also designed to work with both DWTs and MODWTs. universal.thresh <- function (wc, max.level = 4, thres.type = "hard", delta = NULL) { modwt.p <- (class(wc) == "modwt") n <- length(if(modwt.p) wc$d1 else idwt(wc)) wc.fine <- wc[["d1"]] if(is.null(delta)) delta <- median(abs(wc.fine))*sqrt(2*log(n))/0.6745 wc.shrink <- wc thresh.func <- if (thres.type == "hard") function(c) {c * (abs(c) > delta)} else if (thres.type == "soft") function(c) {sign(c) * (abs(c) - delta) * (abs(c) > delta)} else if (thres.type == "mid") function(c) {sign(c) * (abs(c) > delta) * ((abs(c) <= 2*delta) * 2 * (abs(c) - delta) + (abs(c) > 2*delta) * abs(c))} else stop("Invalid option for thres.type (must be \"hard\", \"soft\" or \"mid\")") sqrt2 <- sqrt(2) for (i in names(wc)[1:max.level]) { wc.shrink[[i]] <- thresh.func(wc[[i]]) if(modwt.p) delta <- delta/sqrt2 } wc.shrink } ### read in NMR series using one of the following two R commands nmr <- scan("/Users/dbp/R/Lund/nmr.dat") ### will need to adjust pathname nmr <- scan("http://www.ms.washington.edu/~s530/data/nmr.dat") ### compute the LA(8) DWT nmr.la8.dwt <- dwt(nmr, "la8", 6) ### apply universal threshold with hard, soft and mid thresholding functions nmr.la8.dwt.hard <- universal.thresh(nmr.la8.dwt, 6) nmr.la8.dwt.soft <- universal.thresh(nmr.la8.dwt, 6, thres="soft") nmr.la8.dwt.mid <- universal.thresh(nmr.la8.dwt, 6, thres="mid") ### plot the results for hard thresholding par(mfrow=c(2,1)) plot(nmr,typ="l",ylim=c(-10,60)) enmr.dwt.hard <- idwt(nmr.la8.dwt.hard) plot(enmr.dwt.hard, typ="l", ylim=c(-10,60), main="hard") abline(h=max(nmr), col="red") ### look at effect of thresholding on level 1 ... par(mfrow=c(3,1)) plot(nmr.la8.dwt$d1,typ="h",main="level 1") plot(nmr.la8.dwt.hard$d1,typ="h") plot(nmr.la8.dwt$d1,nmr.la8.dwt.hard$d1) ### on level 2 ... plot(nmr.la8.dwt$d2,typ="h",main="level 2") plot(nmr.la8.dwt.hard$d2,typ="h") plot(nmr.la8.dwt$d2,nmr.la8.dwt.hard$d2) ### on level 3 ... plot(nmr.la8.dwt$d3,typ="h",main="level 3") plot(nmr.la8.dwt.hard$d3,typ="h") plot(nmr.la8.dwt$d3,nmr.la8.dwt.hard$d3) ### on level 4 ... plot(nmr.la8.dwt$d4,typ="h",main="level 4") plot(nmr.la8.dwt.hard$d4,typ="h") plot(nmr.la8.dwt$d4,nmr.la8.dwt.hard$d4) ### on level 5 ... plot(nmr.la8.dwt$d5,typ="h",main="level 5") plot(nmr.la8.dwt.hard$d5,typ="h") plot(nmr.la8.dwt$d5,nmr.la8.dwt.hard$d5) ### and on level 6 plot(nmr.la8.dwt$d6,typ="h",main="level 6") plot(nmr.la8.dwt.hard$d6,typ="h") plot(nmr.la8.dwt$d6,nmr.la8.dwt.hard$d6) ### plot the results for soft thresholding ... par(mfrow=c(2,1)) plot(nmr,typ="l",ylim=c(-10,60)) enmr.dwt.soft <- idwt(nmr.la8.dwt.soft) plot(enmr.dwt.soft, typ="l", ylim=c(-10,60), main="soft") abline(h=max(nmr), col="red") ### ... and for mid thresholding plot(nmr,typ="l",ylim=c(-10,60)) enmr.dwt.mid <- idwt(nmr.la8.dwt.mid) plot(enmr.dwt.mid, typ="l", ylim=c(-10,60), main="mid") abline(h=max(nmr), col="red") ### look at effect of thresholding on level 5 par(mfrow=c(3,1)) plot(nmr.la8.dwt$d5,typ="h",main="level 5") plot(nmr.la8.dwt.mid$d5,typ="h") plot(nmr.la8.dwt$d5,nmr.la8.dwt.mid$d5) ### stack.plot(data.frame(cbind(nmr, enmr.dwt.hard, enmr.dwt.soft, enmr.dwt.mid))) ### compute the LA(8) MODWT nmr.la8.modwt <- modwt(nmr, "la8", 6) ### apply universal threshold with hard, soft and mid thresholding functions nmr.la8.modwt.hard <- universal.thresh(nmr.la8.modwt, 6) nmr.la8.modwt.soft <- universal.thresh(nmr.la8.modwt, 6, thres="soft") nmr.la8.modwt.mid <- universal.thresh(nmr.la8.modwt, 6, thres="mid") ### plot the results for hard thresholding ... par(mfrow=c(2,1)) plot(nmr,typ="l",ylim=c(-10,60)) enmr.modwt.hard <- imodwt(nmr.la8.modwt.hard) plot(enmr.modwt.hard, typ="l", ylim=c(-10,60), main="hard") abline(h=max(nmr), col="red") ### compare with DWT result ... plot(enmr.dwt.hard, typ="l", ylim=c(-10,60), main="hard, DWT") abline(h=max(nmr), col="red") plot(enmr.modwt.hard, typ="l", ylim=c(-10,60), main="hard, MODWT") abline(h=max(nmr), col="red") ### ... for soft thresholding ... plot(nmr,typ="l",ylim=c(-10,60)) enmr.modwt.soft <- imodwt(nmr.la8.modwt.soft) plot(enmr.modwt.soft, typ="l", ylim=c(-10,60), main="soft") abline(h=max(nmr), col="red") ### compare with DWT result ... plot(enmr.dwt.soft, typ="l", ylim=c(-10,60), main="soft, DWT") abline(h=max(nmr), col="red") plot(enmr.modwt.soft, typ="l", ylim=c(-10,60), main="soft, MODWT") abline(h=max(nmr), col="red") ### ... and for mid thresholding plot(nmr,typ="l",ylim=c(-10,60)) enmr.modwt.mid <- imodwt(nmr.la8.modwt.mid) plot(enmr.modwt.mid, typ="l", ylim=c(-10,60), main="mid") abline(h=max(nmr), col="red") ### compare with DWT result ... plot(enmr.dwt.mid, typ="l", ylim=c(-10,60), main="mid, DWT") abline(h=max(nmr), col="red") plot(enmr.modwt.mid, typ="l", ylim=c(-10,60), main="mid, MODWT") abline(h=max(nmr), col="red") ### stack.plot(data.frame(cbind(nmr, enmr.modwt.hard, enmr.modwt.soft, enmr.modwt.mid))) ### read in ECG series using one of the following two R commands ecg <- scan("/Users/dbp/R/Lund/heart.dat") ### will need to adjust pathname ecg <- scan("http://www.ms.washington.edu/~s530/data/heart.dat") ### compute the LA(8) DWT ecg.la8.dwt <- dwt(ecg, "la8", 6) ### apply universal threshold with hard, soft and mid thresholding functions ### and set scaling coefficients to zero rather than keeping them unaltered ecg.la8.dwt.hard <- universal.thresh(ecg.la8.dwt, 6) ecg.la8.dwt.hard$s6 <- rep(0,length(ecg.la8.dwt.hard$s6)) ecg.la8.dwt.soft <- universal.thresh(ecg.la8.dwt, 6, thres="soft") ecg.la8.dwt.soft$s6 <- rep(0,length(ecg.la8.dwt.soft$s6)) ecg.la8.dwt.mid <- universal.thresh(ecg.la8.dwt, 6, thres="mid") ecg.la8.dwt.mid$s6 <- rep(0,length(ecg.la8.dwt.mid$s6)) ### plot the results for hard thresholding par(mfrow=c(2,1)) plot(ecg,typ="l",ylim=c(-2,2)) eecg.dwt.hard <- idwt(ecg.la8.dwt.hard) plot(eecg.dwt.hard, typ="l", ylim=c(-2,2), main="hard") ### look at effect of thresholding on level 1 ... par(mfrow=c(3,1)) plot(ecg.la8.dwt$d1,typ="h",main="level 1") plot(ecg.la8.dwt.hard$d1,typ="h") plot(ecg.la8.dwt$d1,ecg.la8.dwt.hard$d1) ### on level 2 ... plot(ecg.la8.dwt$d2,typ="h",main="level 2") plot(ecg.la8.dwt.hard$d2,typ="h") plot(ecg.la8.dwt$d2,ecg.la8.dwt.hard$d2) ### on level 3 ... plot(ecg.la8.dwt$d3,typ="h",main="level 3") plot(ecg.la8.dwt.hard$d3,typ="h") plot(ecg.la8.dwt$d3,ecg.la8.dwt.hard$d3) ### on level 4 ... plot(ecg.la8.dwt$d4,typ="h",main="level 4") plot(ecg.la8.dwt.hard$d4,typ="h") plot(ecg.la8.dwt$d4,ecg.la8.dwt.hard$d4) ### on level 5 ... plot(ecg.la8.dwt$d5,typ="h",main="level 5") plot(ecg.la8.dwt.hard$d5,typ="h") plot(ecg.la8.dwt$d5,ecg.la8.dwt.hard$d5) ### and on level 6 plot(ecg.la8.dwt$d6,typ="h",main="level 6") plot(ecg.la8.dwt.hard$d6,typ="h") plot(ecg.la8.dwt$d6,ecg.la8.dwt.hard$d6) ### plot the results for soft thresholding ... par(mfrow=c(2,1)) plot(ecg,typ="l",ylim=c(-2,2)) eecg.dwt.soft <- idwt(ecg.la8.dwt.soft) plot(eecg.dwt.soft, typ="l", ylim=c(-2,2), main="soft") ### compare hard and soft plot(eecg.dwt.hard, typ="l", ylim=c(-2,2), main="hard") abline(h=1.0, col="red") abline(h=-1.0, col="red") plot(eecg.dwt.soft, typ="l", ylim=c(-2,2), main="soft") abline(h=1.0, col="red") abline(h=-1.0, col="red") ### ... and for mid thresholding plot(ecg,typ="l",ylim=c(-2,2)) eecg.dwt.mid <- idwt(ecg.la8.dwt.mid) plot(eecg.dwt.mid, typ="l", ylim=c(-2,2), main="mid") ### compare hard and mid plot(eecg.dwt.hard, typ="l", ylim=c(-2,2), main="hard") abline(h=1.0, col="red") abline(h=-1.0, col="red") plot(eecg.dwt.mid, typ="l", ylim=c(-2,2), main="mid") abline(h=1.0, col="red") abline(h=-1.0, col="red") ### stack.plot(data.frame(cbind(ecg, eecg.dwt.hard, eecg.dwt.soft, eecg.dwt.mid))) ### compute the LA(8) MODWT ecg.la8.modwt <- modwt(ecg, "la8", 6) ### apply universal threshold with hard, soft and mid thresholding functions ecg.la8.modwt.hard <- universal.thresh(ecg.la8.modwt, 6) ecg.la8.modwt.hard$s6 <- rep(0,length(ecg.la8.modwt.hard$s6)) ecg.la8.modwt.soft <- universal.thresh(ecg.la8.modwt, 6, thres="soft") ecg.la8.modwt.soft$s6 <- rep(0,length(ecg.la8.modwt.soft$s6)) ecg.la8.modwt.mid <- universal.thresh(ecg.la8.modwt, 6, thres="mid") ecg.la8.modwt.mid$s6 <- rep(0,length(ecg.la8.modwt.mid$s6)) ### plot the results for hard thresholding ... par(mfrow=c(2,1)) plot(ecg,typ="l",ylim=c(-2,2)) eecg.modwt.hard <- imodwt(ecg.la8.modwt.hard) plot(eecg.modwt.hard, typ="l", ylim=c(-2,2), main="hard") ### compare with DWT result ... plot(eecg.dwt.hard, typ="l", ylim=c(-2,2), main="hard, DWT") abline(h=1.0, col="red") abline(h=-1.0, col="red") plot(eecg.modwt.hard, typ="l", ylim=c(-2,2), main="hard, MODWT") abline(h=1.0, col="red") abline(h=-1.0, col="red") ### ... for soft thresholding ... plot(ecg,typ="l",ylim=c(-2,2)) eecg.modwt.soft <- imodwt(ecg.la8.modwt.soft) plot(eecg.modwt.soft, typ="l", ylim=c(-2,2), main="soft") ### compare with DWT result ... plot(eecg.dwt.soft, typ="l", ylim=c(-2,2), main="soft, DWT") abline(h=1.0, col="red") abline(h=-1.0, col="red") plot(eecg.modwt.soft, typ="l", ylim=c(-2,2), main="soft, MODWT") abline(h=1.0, col="red") abline(h=-1.0, col="red") ### ... and for mid thresholding plot(ecg,typ="l",ylim=c(-2,2)) eecg.modwt.mid <- imodwt(ecg.la8.modwt.mid) plot(eecg.modwt.mid, typ="l", ylim=c(-2,2), main="mid") ### compare with DWT result ... plot(eecg.dwt.mid, typ="l", ylim=c(-2,2), main="mid, DWT") abline(h=1.0, col="red") abline(h=-1.0, col="red") plot(eecg.modwt.mid, typ="l", ylim=c(-2,2), main="mid, MODWT") abline(h=1.0, col="red") abline(h=-1.0, col="red") ### stack.plot(data.frame(cbind(ecg, eecg.modwt.hard, eecg.modwt.soft, eecg.modwt.mid)))