threshsure <- function(ywd, verbose = F)
{
  ans <- ywd
  n <- length(ywd$D)
  nlev <- log(n + 1, base = 2) - 1
  i <- nlev
  iloc <- 1
  while(i >= 0) {
#       Extract current level coefficients from all wavelet coefficients
    gg <- sort(abs(ywd$D[iloc:(iloc + 2^i - 1)]))	
#       Calculate SURE for all possible thresholds
    surevec <- sapply(gg, sure, gg)	
#       Threshold that minimizes risk
    if(sure(0, gg) < min(surevec))
    else
	thresh <- gg[match(min(surevec), surevec)]
    if(verbose)
      cat(paste("At level ", i, ", the threshold is ", thresh,
		"\n", sep = ""))
#       Perform shrinking
    ans$D[iloc:(iloc + 2^i - 1)] <- shrinkit(ywd$D[iloc:(iloc + 2^i -
							 1)], thresh)
    iloc <- iloc + 2^i
    i <- i - 1
  }
  ans
}

shrinkit <- function(coeffs, thresh)
{
  sign(coeffs) * pmax(abs(coeffs) - thresh, 0)
}

sure <- function(t, x)
{
  ax <- sort(abs(x))	
					#  num gives the number of elements
					# in x that are <= t
  num <- match(F, ax <= t, nomatch = length(ax) + 1) - 1
  length(ax) - 2 * num + sum(pmin(ax, t)^2)
}