R/denoiseheteroprop.R

In adlift: An Adaptive Lifting Scheme Algorithm

"denoiseheteroprop" <-function(x,f,pred=LinearPred,neigh=1,int=TRUE,clo=FALSE,keep=2,rule="median",gamvec=rep(1,times=length(x)),returnall=FALSE){

#denoises input using prediction scheme specified by other arguments
#artificial level thresholding, values generated by ebayesthresh using normalized detail coefficients

#with info about variances:gamvec is a vector of noise sd (!SD, not var!) proportionality constants
# (in order of observing x)

newcoeff<-NULL
ndetlist<-list()
tclist<-NULL
norcoeff<-NULL

n<-length(x)

out<-fwtnp(x,f,LocalPred=pred,neighbours=neigh,intercept=int,closest=clo,nkeep=keep,do.W=TRUE,varonly=FALSE)
w<-out$W


#work out sd for normalising...

sigmamat<-diag(gamvec^2)	     #sigma is a diag matrix with sigmai^2=gamvec[i] on the leading diag

varmat<-w%*%(sigmamat)	     
varmat<-(varmat)%*%t(w)	     #varmat is the var-covariance matrix of the lifted coeff's
indsd<-sqrt(diag(varmat))    #this will contain the sd of the scaling coefficients as well

norcoeff<-out$coeff/indsd

nonorcoeff<-out$coeff

lr<-out$lengthsremove    #vector deciding how to divide up coefficients into artificial levels
rem<-out$removelist      #used to convert output to original lr,rem)
po<-out$pointsin

al<-artlev(lr,rem)      #the list of indices of removelist separated into levels

levno<-length(al)
for (i in 1:levno){
	ndetlist[[i]]<-norcoeff[al[[i]]]
}

sd<-mad(ndetlist[[1]])    #uses the first (largest) level to estimate noise standard deviation
sd1<-mad(ndetlist[[1]],center=0)
sd2<-mad(norcoeff[rem]) #uses all the normalised details to estimate the noise std dev

for (i in 1:levno){
	tclist<-ebayesthresh(ndetlist[[i]],prior="cauchy",a=NA,sdev=sd,threshrule=rule)
	newcoeff[al[[i]]]<-tclist*indsd[al[[i]]]	#creates list of thresholded coeffs (levels) using HARD thresholding and post. median values 
}

newcoeff[po]<-out$coeff[po]

fhat<-invtnp(x,newcoeff,out$lengths,lr,po,rem,out$neighbrs,out$schemehist,out$interhist,length(x)-keep,int,neigh,clo,pred)

if(returnall){
	return(list(fhat=fhat,w=w,indsd=indsd,al=al,sd=sd))
}
else{
	return(fhat$coeff)
}
}