bitup.ti <- function(noisy, filter.number, family) {
	n <- length(noisy)
	tmp1 <- parent.ti(noisy, filter.number, family)
	tmp2 <- parent.ti(noisy[n:1], filter.number, family)
	tmp2 <- tmp2[n:1]
	return((tmp1 + tmp2)/2)
}

bitup.ti.shift <- function(noisy, filter.number, family, max.shift = 3) {
	n <- length(noisy)
	tmp1 <- parent.ti.shift(noisy, filter.number, family, max.shift)
	tmp2 <- parent.ti.shift(noisy[n:1], filter.number, family, max.shift)
	tmp2 <- tmp2[n:1]
	return((tmp1 + tmp2)/2)
}

shift.sequence <- function(vect, m)
{
	len <- length(vect)
	if(m >= 1.)
		return(c(vect[(m + 1.):len], vect[1.:m]))
	else if(m <= -1.)
		return(c(vect[(len + m + 1.):len], vect[1.:(len + m)]))
	else return(vect)
}

scale.indices <- function(N, scale, type = "decimated")
{
	if(type == "decimated") {
		beg <- 1
		endd <- N/2
		if(scale > 1)
			for(i in 1:(scale - 1)) {
				beg.old <- beg
				beg <- endd + 1
				endd <- (endd - beg.old + 1)/2 + endd
			}
	}
	if(type == "station") {
		beg <- (scale - 1) * N + 1
		endd <- scale * N
	}
	return(beg, endd)
}

which.shift <- function(a, b, max.jumps, length.jump, power = 1) {
	n <- length(a)
	max.jumps <- min(max.jumps, floor((n-1)/length.jump))
	al <- rep(0, 2*max.jumps+1)
	for (i in (-max.jumps):max.jumps) al[i+max.jumps+1] <- sum(abs(abs(a) - abs(shift.sequence(b, i*length.jump)))^power)
	return((which(al == min(al)) - max.jumps - 1)*length.jump)
}

parent.ti <- function(noisy, filter.number, family) {
	n <- length(noisy)
	J <- logb(n, 2)
	noisy.w <- wd(noisy, filter.number, family, type="station")
	sigma <- (mad(noisy.w$D[2 * (1:(n/2))]) + mad(noisy.w$D[2 * (1:(n/2)) - 1]))/2
	tmp <- noisy.w
	for (i in 1:(J-1)) {
		v <- scale.indices(n, i, "station")
		v1 <- scale.indices(n, i+1, "station")
		indices <- rep((2^i+1):(2^(i+1)), n * 2^(-i-1)) + rep(seq(from=0, by = 2^(i+1), length=n * 2^(-i-1)), each = 2^i)
		temp <- tmp$D[v1$beg:v1$endd]
		temp[indices] <- temp[indices - 2^i]
		tmp$D[v$beg:v$endd] <- tmp$D[v$beg:v$endd]^2 + temp^2
	}
	v <- scale.indices(n, J, "station")
	tmp$D[v$beg:v$endd] <- tmp$D[v$beg:v$endd]^2
	noisy.w$D <- noisy.w$D * (tmp$D > sigma^2 * 2 * log(n))
	return(AvBasis(convert((noisy.w))))
}

parent.ti.shift <- function(noisy, filter.number, family, max.shift = 3) {
	n <- length(noisy)
	J <- logb(n, 2)
	noisy.w <- wd(noisy, filter.number, family, type="station")
	sigma <- (mad(noisy.w$D[2 * (1:(n/2))]) + mad(noisy.w$D[2 * (1:(n/2)) - 1]))/2
	tmp <- noisy.w
	for (i in 1:(J-1)) {
		v <- scale.indices(n, i, "station")
		v1 <- scale.indices(n, i+1, "station")
		indices <- rep((2^i+1):(2^(i+1)), n * 2^(-i-1)) + rep(seq(from=0, by = 2^(i+1), length=n * 2^(-i-1)), each = 2^i)
		temp <- tmp$D[v1$beg:v1$endd]
		temp[indices] <- temp[indices - 2^i]
		ss <- median(which.shift(tmp$D[v$beg:v$endd], temp, max.shift, 2^(i+1)))
		tmp$D[v$beg:v$endd] <- tmp$D[v$beg:v$endd]^2 + shift.sequence(temp, ss)^2
	}
	v <- scale.indices(n, J, "station")
	tmp$D[v$beg:v$endd] <- tmp$D[v$beg:v$endd]^2
	noisy.w$D <- noisy.w$D * (tmp$D > sigma^2 * 2 * log(n))
	return(AvBasis(convert((noisy.w))))
}

######################################################

babte.ti <- function(noisy, w1, w2 = w1, family1 = "DaubExPhase", family2 = "DaubLeAsymm") {

	TT <- length(noisy)
	J <- logb(TT, 2)

	noisy.w1 <- wd(noisy, w1, family1, type="station")
	noisy.w2 <- wd(noisy, w2, family2, type="station")

	sigma <- (mad(noisy.w1$D[2*(1:(TT/2))]) + mad(noisy.w1$D[2*(1:(TT/2))-1]) + mad(noisy.w2$D[2*(1:(TT/2))]) + 
mad(noisy.w2$D[2*(1:(TT/2))-1]))/4
	
	weights <- rep(J*TT, 0)

	f1 <- hp(w1, family1, J)
	f2 <- hp(w2, family2, J)
	
	for (j in 1:J) {
	
		A <- matrix(sigma^2, 2, 2)
		overlap.length <- min(length(f1[[j]]), length(f2[[j]]))
		A[2,1] <- A[1,2] <- sigma^2 * sum(f1[[j]][(1:overlap.length)+2^j*max(w1-w2,0)] * 
f2[[j]][(1:overlap.length)+2^j*max(w2-w1,0)])

		B <- solve(A)

		indices <- scale.indices(TT, j, "station")
		
		D <- matrix(0, TT, 2)
		D[,1] <- noisy.w1$D[indices$beg:indices$endd]
		D[,2] <- noisy.w2$D[indices$beg:indices$endd]

		E <- (D %*% B) * D
	
		FF <- (E[,1] + E[,2])

		weights[indices$beg:indices$endd] <- (FF > 2 * log(TT))

	}

	rec.joint.w1 <- noisy.w1
	rec.joint.w2 <- noisy.w2
	
	rec.joint.w1$D <- rec.joint.w1$D * weights
	rec.joint.w2$D <- rec.joint.w2$D * weights
	
	rec.joint <- (AvBasis(convert(rec.joint.w1)) + AvBasis(convert(rec.joint.w2)))/2
	
	return(rec.joint)
	
}

gname <- function(filter.number, filter.family, max.scale = 1) {
	return(paste("g.",max.scale,".",filter.number,".",filter.family, sep=""))
}

hp <- function(filter.number, filter.family, max.scale = 1) {
	gn <- gname(filter.number, filter.family, max.scale)
	if (exists(gn)) return(get(gn))

	cat("Creating discrete wavelet vectors, please wait...\n")

	res <- vector("list", max.scale)
	
	tmp <- highpass(filter.number, filter.family)
	res[[1]] <- tmp
	h <- filter.select(filter.number, filter.family)$H

	if (max.scale > 1) for (i in 2:max.scale) {
		wavlength.old <- length(res[[i-1]])
		wavlength <- (2^i - 1) * (filter.number*2 - 1) + 1
		tmp <- rep(0, wavlength)
		
		for (n in 0:(wavlength-1)) {
			k <- (max(0, ceiling(n/2 - filter.number + 1/2))) : (min(wavlength.old-1, floor(n/2)))
			tmp[n+1] <- sum(h[n-2*k+1] * res[[i-1]][k+1])
		}
		res[[i]] <- tmp
	}
	assign(gn, res, where = 1)

	cat("...done.\n")

	return(res)
}

highpass <- function(filter.number, filter.family) {
	res <- filter.select(filter.number, filter.family)$H
	N <- length(res)
	res <- res[N:1] * (-1)^(0:(N-1))
	return(res)
}