sqrt.est <- function(noisy, wav.filter.number = 4, family = "DaubLeAsymm", theta = 0.01, ti = T) {
	n <- length(noisy)
	J <- logb(n, 2)

	if (ti) {
		noisy.w <- wd(noisy, wav.filter.number, family, type="station")
		sigma <- (mad(noisy.w$D[2*(1:(n/2))]) + mad(noisy.w$D[2*(1:(n/2))-1]))/2
	}
	else {
		noisy.w <- wd(noisy, wav.filter.number, family)
		sigma <- mad(noisy.w$D[1:(n/2)])
	}

	thresholds <- sqrt(theta + (1-theta)*(0:(J-1))/(J-1)) * sigma * sqrt(J * 2 * log(2))
	
	noisy.w.t <- threshold(noisy.w, policy="manual", type="hard", value = thresholds, levels=0:(noisy.w$nlevels-1))

	if (ti)
		noisy.w.t.r <- AvBasis(convert(noisy.w.t))
	else
		noisy.w.t.r <- wr(noisy.w.t)
		
	return(noisy.w.t.r)
}

cv.score.theta <- function(noisy, theta, wav.filter.number = 1, family = "DaubExPhase") {
	n <- length(noisy)
	odd.indices <- seq(from = 1, by = 2, to = n)
	even.indices <- seq(from = 2, by = 2, to = n)
	noisy.odd <- noisy[odd.indices]
	noisy.even <- noisy[even.indices]
	noisy.odd.rec <- new.univ(noisy.odd, wav.filter.number, family, theta, ti = F)
	noisy.even.rec <- new.univ(noisy.even, wav.filter.number, family, theta, ti = F)
	return(sum((noisy.even.rec - noisy.odd)^2) + sum((noisy.odd.rec - noisy.even)^2))
}

find.theta <- function(noisy, thetamin = 0, thetamax = 1, thetalength = 11, wav.filter.number = 1, family = "DaubExPhase") {
	cvsc <- rep(0, thetalength)
	theta <- seq(from = thetamin, to = thetamax, length = thetalength)
	for (i in 1:thetalength) cvsc[i] <- cv.score.theta(noisy, theta[i], wav.filter.number, family)
	min.index <- which(cvsc == min(cvsc))
	mil <- length(min.index)
	if (min.index[mil] == thetalength) min.index <- thetalength else
		if (min.index[1] == 1) min.index <- 1
			else min.index <- round(median(min.index))
	return(theta[min.index])
}

sqrt.est.cv <- function(noisy, wav.filter.number = 4, family = "DaubLeAsymm", thetamin = 0.01, thetamax = 1, thetalength = 11, ti = T) {
	theta <- find.theta(noisy, thetamin, thetamax, thetalength, wav.filter.number, family)
	return(sqrt.est(noisy, wav.filter.number, family, theta, ti))
}