Move baselineWavelet to github from google code

DESCRIPTION

@@ -0,0 +1,12 @@
+Package: baselineWavelet
+Type: Package
+Title: baseline correction by wavelet and Whittaker Smooth algorithms
+Version: 4.0.1
+Date: 2012-11-23
+Depends: R (>= 2.15.0), Matrix
+Suggests:
+Author: Zhimin Zhang, Yizeng Liang
+Maintainer: Zhimin Zhang <[email protected]>
+Description: baseline correction by wavelet and Whittaker Smooth algorithms
+License: LGPL version 2 or newer
+Packaged: Fri Nov 23 20:40:21 2012;

R/WhittakerSmooth.R

@@ -0,0 +1,9 @@
+WhittakerSmooth <- function(x,w,lambda,differences=1) {
+  x=as.vector(x)
+  L=length(x)
+  E=spMatrix(L,L,i=seq(1,L),j=seq(1,L),rep(1,L))
+  D=as(diff(E,1,differences),"dgCMatrix")
+  W=as(spMatrix(L,L,i=seq(1,L),j=seq(1,L),w),"dgCMatrix")
+  background=solve((W+lambda*t(D)%*%D),w*x);
+  return(as.vector(background))
+ }

R/baselineCorrectionCWT.R

@@ -0,0 +1,114 @@
+"baselineCorrectionCWT" <-
+function(x,peakWidth,threshold=0.5,lambda=100,differences=1) {
+
+  # fitting an rough background with proper value of lambda
+  x.w=rep(1,length(x))
+  peakIndex=peakWidth$peakIndex
+  for (i in 1:length(peakIndex)){
+    x.w[peakWidth[[paste(peakIndex[i])]]]=0
+  }
+  backgr = WhittakerSmooth(x,x.w,lambda,differences)
+
+  #assigning the corresponding value of the spectrum to the part of the rough background,which is larger than the spectrum. and fitting another background without value lager than orignal spectrum based on the new rough background
+
+  backgr.final=backgr
+  backgr.final[x<=backgr]=x[x<=backgr]
+  backgr.final= WhittakerSmooth(backgr.final,rep(1,length(backgr.final)),differences)
+  x= backgr.final
+
+  #refine the new rough background
+
+  background = NULL
+  peakIndex=peakWidth$peakIndex
+  signal_baseline= x[1:(peakWidth[[paste(peakIndex[1])]][1])-1]
+  signal_baseline.w= rep(1,length(signal_baseline))
+  baseline_baseline =   WhittakerSmooth(signal_baseline,signal_baseline.w,differences)
+  background = c(background,baseline_baseline)
+  for (i in 1:(length(peakIndex)-1)){
+
+    peakWidth.1=peakWidth[[paste(peakIndex[i])]]
+    peakWidth.2=peakWidth[[paste(peakIndex[i+1])]]
+
+    if(length(intersect(peakWidth.1, peakWidth.2))==0)
+    {
+      signal_peak=x[peakWidth.1]
+      signal_peak.w=c(1,rep(0,(length(peakWidth.1)-2)),1)
+      if((abs(signal_peak[length(signal_peak)]-signal_peak[1])/mean(signal_peak-min(signal_peak)))>=threshold){
+        if(signal_peak[length(signal_peak)]>signal_peak[1]){
+          signal_peak.w=as.numeric(!(signal_peak>signal_peak[length(signal_peak)]))
+        }else{
+          signal_peak.w=as.numeric(!(signal_peak>signal_peak[1]))
+        }
+        signal_peak.w[1]=1
+        signal_peak.w[length(signal_peak.w)]=1
+      }
+      peak_baseline= WhittakerSmooth(signal_peak,signal_peak.w,differences)
+      background = c(background,peak_baseline)
+
+      if(peakWidth.2[1]-peakWidth.1[length(peakWidth.1)]==1){
+        # two peak just together
+      }else{
+        signal_baseline=x[(peakWidth.1[length(peakWidth.1)]+1):(peakWidth.2[1]-1)]
+        if(length(signal_baseline)<=3){
+          baseline_baseline=signal_baseline
+        }else{
+          signal_baseline.w=c(1,rep(1,(length(signal_baseline)-2)),1)
+	    baseline_baseline= WhittakerSmooth(signal_baseline,signal_baseline.w,differences)
+        }
+        background = c(background,baseline_baseline)
+      }
+
+    }else{
+      if(peakWidth.1[1]>peakWidth.2[1]){
+         # the last peak contain the preivous peak
+         peakWidth[[paste(peakIndex[i+1])]]=(peakWidth.1[1]):(peakWidth.2[length(peakWidth.2)])
+      }else{
+        if(length(peakWidth.1)>=length(peakWidth.2)){
+          peakWidth.11=setdiff(peakWidth.1,intersect(peakWidth.1, peakWidth.2))
+        }else{
+          peakWidth.11=peakWidth.1
+          peakWidth[[paste(peakIndex[i+1])]]=setdiff(peakWidth.2,intersect(peakWidth.1, peakWidth.2))
+        }
+
+        if(length(peakWidth.11)<=2){
+          peakWidth[[paste(peakIndex[i+1])]]=(peakWidth.1[1]):(peakWidth.2[length(peakWidth.2)])
+        } else{
+          signal_peak=x[peakWidth.11]
+          signal_peak.w=c(1,rep(0,(length(peakWidth.11)-2)),1)
+          if((abs(signal_peak[length(signal_peak)]-signal_peak[1])/mean(signal_peak-min(signal_peak)))>=threshold){
+            if(signal_peak[length(signal_peak)]>signal_peak[1]){
+              signal_peak.w=as.numeric(!(signal_peak>signal_peak[length(signal_peak)]))
+            }else{
+              signal_peak.w=as.numeric(!(signal_peak>signal_peak[1]))
+            }
+            signal_peak.w[1]=1
+            signal_peak.w[length(signal_peak.w)]=1
+          }
+          peak_baseline= WhittakerSmooth(signal_peak,signal_peak.w,differences)
+          background = c(background,peak_baseline)
+        }
+
+      }
+    }
+  }
+
+  peakWidth.end= peakWidth[[paste(peakIndex[length(peakIndex)])]]
+  signal_peak=x[peakWidth.end]
+  signal_peak.w=c(1,rep(0,(length(peakWidth.end)-2)),1)
+  if((abs(signal_peak[length(signal_peak)]-signal_peak[1])/mean(signal_peak-min(signal_peak)))>=threshold){
+    if(signal_peak[length(signal_peak)]>signal_peak[1]){
+      signal_peak.w=as.numeric(!(signal_peak>signal_peak[length(signal_peak)]))
+    }else{
+      signal_peak.w=as.numeric(!(signal_peak>signal_peak[1]))
+    }
+  }
+  peak_baseline= WhittakerSmooth(signal_peak,signal_peak.w,1)
+
+  signal_baseline=x[(peakWidth.end[length(peakWidth.end)]+1):length(x)]
+  signal_baseline.w=rep(1,length(signal_baseline))
+  baseline_baseline= WhittakerSmooth(signal_baseline,signal_baseline.w,1)
+  background = c(background,peak_baseline)
+  background = c(background,baseline_baseline)
+
+  return(background)
+}

R/cwt.R

@@ -0,0 +1,55 @@
+"cwt" <-
+function(ms, scales=1, wavelet='mexh') {
+	## Check for the wavelet format
+	if (wavelet == 'mexh') {
+		psi_xval <- seq(-8, 8, length=1024)
+		psi <- (2/sqrt(3) * pi^(-0.25)) * (1 - psi_xval^2) *exp(-psi_xval^2/2)
+		#plot(psi_xval, psi)
+	} else if (wavelet=='haar') {
+    psi_xval <- seq(0,1,length=1024)
+    psi <- c(0,rep(1,511),rep(-1,511),0)
+	}else if (is.matrix(wavelet)) {
+		if (nrow(wavelet) == 2) {
+			psi_xval <- wavelet[1,]
+			psi <- wavelet[2,]
+		} else if (ncol(wavelet) == 2) {
+			psi_xval <- wavelet[,1]
+			psi <- wavelet[,2]
+		} else {
+			stop('Unsupported wavelet format!')
+		}
+	} else {
+		stop('Unsupported wavelet!')
+	}
+
+    oldLen <- length(ms)
+	## To increase the computation effeciency of FFT, extend it as the power of 2
+	## because of a strange signal length 21577 makes the FFT very slow!
+	#ms <- extend.nBase(ms, nLevel=1, base=2)
+	ms <- extendNBase(ms, nLevel=NULL, base=2)
+	len <- length(ms)
+    nbscales <- length(scales)
+    wCoefs <- NULL
+
+    psi_xval <- psi_xval - psi_xval[1]
+    dxval <- psi_xval[2]
+    xmax  <- psi_xval[length(psi_xval)]
+    for (i in 1:length(scales)) {
+		scale.i <- scales[i]
+		f <- rep(0, len)
+        j <- 1 + floor((0:(scale.i * xmax))/(scale.i * dxval))
+        if (length(j) == 1)		j <- c(1, 1)
+		lenWave <- length(j)
+        f[1:lenWave] <- rev(psi[j]) - mean(psi[j])
+		if (length(f) > len) stop(paste('scale', scale.i, 'is too large!'))
+		wCoefs.i <- 1/sqrt(scale.i) * convolve(ms, f)
+		## Shift the position with half wavelet width
+		wCoefs.i <- c(wCoefs.i[(len-floor(lenWave/2) + 1) : len], wCoefs.i[1:(len-floor(lenWave/2))])
+		wCoefs <- cbind(wCoefs, wCoefs.i)
+    }
+	if (length(scales) == 1) wCoefs <- matrix(wCoefs, ncol=1)
+	colnames(wCoefs) <- scales
+	wCoefs <- wCoefs[1:oldLen,,drop=FALSE]
+	return(wCoefs)
+}
+

R/extendLength.R

@@ -0,0 +1,37 @@
+"extendLength" <-
+function(x, addLength=NULL, method=c('reflection', 'open', 'circular'), direction=c('right', 'left', 'both')) {
+	if (is.null(addLength)) stop('Please provide the length to be added!')
+	if (!is.matrix(x)) x <- matrix(x, ncol=1)	
+	method <- match.arg(method)
+	direction <- match.arg(direction)
+
+	nR <- nrow(x)
+	nR1 <- nR + addLength
+	if (direction == 'both') {
+		left <- right <- addLength
+	} else if (direction == 'right') {
+		left <- 0
+		right <- addLength
+	} else if (direction == 'left') {
+		left <- addLength
+		right <- 0
+	}
+
+	if (right > 0) {
+		x <- switch(method,
+			reflection =rbind(x, x[nR:(2 * nR - nR1 + 1), , drop=FALSE]),
+			open = rbind(x, matrix(rep(x[nR,], addLength), ncol=ncol(x), byrow=TRUE)),
+			circular = rbind(x, x[1:(nR1 - nR),, drop=FALSE]))
+	}
+
+	if (left > 0) {
+		x <- switch(method,
+			reflection =rbind(x[addLength:1, , drop=FALSE], x),
+			open = rbind(matrix(rep(x[1,], addLength), ncol=ncol(x), byrow=TRUE), x),
+			circular = rbind(x[(2 * nR - nR1 + 1):nR,, drop=FALSE], x))
+	}
+	if (ncol(x) == 1)  x <- as.vector(x)
+
+	return(x)
+}
+

R/extendNBase.R

@@ -0,0 +1,17 @@
+"extendNBase" <-
+function(x, nLevel=1, base=2, ...) {
+	if (!is.matrix(x)) x <- matrix(x, ncol=1)	
+
+	nR <- nrow(x)
+	if (is.null(nLevel)) {
+		nR1 <- nextn(nR, base)		
+	} else {
+		nR1 <- ceiling(nR / base^nLevel) * base^nLevel		
+	}
+	if (nR != nR1) {
+		x <- extendLength(x, addLength=nR1-nR, ...)
+	}
+
+	return(x)
+}
+

R/getLocalMaximumCWT.R

@@ -0,0 +1,21 @@
+"getLocalMaximumCWT" <-
+function(wCoefs, minWinSize=5, amp.Th=0) {
+
+	localMax <- NULL
+	scales <- as.numeric(colnames(wCoefs))
+
+	for (i in 1:length(scales)) {
+		scale.i <- scales[i]
+		winSize.i <- scale.i * 2 + 1
+		if (winSize.i < minWinSize) {
+			winSize.i <- minWinSize
+		} 
+		temp <- localMaximum(wCoefs[,i], winSize.i)
+		localMax <- cbind(localMax, temp)
+	}
+	# Set the values less than peak threshold as 0
+	localMax[wCoefs < amp.Th] <- 0
+	colnames(localMax) <- colnames(wCoefs)
+	rownames(localMax) <- rownames(wCoefs)
+	return(localMax)
+}