filter DO better for nighttime regression

[aappling-usgs]

rather than using filter(), use a sevitsky golay filter (more appropriate to curved time series; then could use the derivatives directly from the polynomial filter). their paper has been cited lots: Smoothing and Differentiation of Data by Simplified Least Squares Procedures. sgolay()

Code from Bob, 9/3/15:

library(signal)

nightregsg<-function(o2file, bp, ts){

  temp<-o2file$temp
  oxy<-o2file$oxy

  ##filter oxy data.  No derivative!
  oxyf<-sgolayfilt(x=oxy, p=3, n=5, m=0)

  ##calculate delO/delt, SG filter with m=1 so we took derivative
  deltaO2<-1440*(sgolayfilt(x=oxy, p=3, n=5, m=1))/ts

  #calc the dodef
  satdef<-osat(temp,bp)-oxyf

  #calculate regression and plot
  nreg<-lm(deltaO2~satdef)
  plot(satdef,deltaO2)
  abline(nreg)

  coeff<-coef(nreg)

  out<-list(coeff, K600fromO2(mean(temp), coeff[2]))
  out
}

nightregsg(french[french$dtime>=as.numeric(chron(dates="09/15/12", times="19:40:00")) & french$dtime<=as.numeric(chron(dates="09/15/12", times="23:00:00")), ], bp=523, ts=5)

nightregsg(french[french$dtime>=as.numeric(chron(dates="08/24/12", times="19:40:00")) & french$dtime<=as.numeric(chron(dates="08/24/12", times="23:00:00")), ], bp=523, ts=5)

nightregsg(french[french$dtime>=as.numeric(chron(dates="08/25/12", times="19:40:00")) & french$dtime<=as.numeric(chron(dates="08/25/12", times="23:00:00")), ], bp=523, ts=5)


##ignore the below  

o2file<-french[french$dtime>=as.numeric(chron(dates="09/15/12", times="19:40:00")) & french$dtime<=as.numeric(chron(dates="09/15/12", times="23:00:00")), ] 
plot(french15$oxy)
length(french15$oxy)

plot(sgolayfilt(x=french15$oxy, p=3, n=7, m=0))
plot(sgolayfilt(x=french15$oxy, p=3, n=7, m=1))
length(sgolayfilt(x=french15$oxy, p=3, n=7, m=1))

[aappling-usgs]

More code from bob, 9/3/15:

# ---------------------------------------------------------------------- 
# Savitzky-Golay Algorithm 
# ---------------------------------------------------------------------- 
# T2 <- sav.gol(T, fl, forder=4, dorder=0); 
# 
# Polynomial filtering method of Savitzky and Golay 
# See Numerical Recipes, 1992, Chapter 14.8, for details. 
# 
# T = vector of signals to be filtered 
# (the derivative is calculated for each ROW) 
# fl = filter length (for instance fl = 51..151) 
# forder = filter order (2 = quadratic filter, 4= quartic) 
# dorder = derivative order (0 = smoothing, 1 = first derivative, etc.) 
# 
sav.gol <- function(T, fl, forder=4, dorder=0) 
{ 
    m <- length(T) 
    dorder <- dorder + 1 
    # -- calculate filter coefficients -- 
    fc <- (fl-1)/2 # index: window left and right 
    X <- outer(-fc:fc, 0:forder, FUN="^") # polynomial terms and coefficients 
    Y <- pinv(X); # pseudoinverse 
    # -- filter via convolution and take care of the end points -- 
    T2 <- convolve(T, rev(Y[dorder,]), type="o") # convolve(...) 
    T2 <- T2[(fc+1):(length(T2)-fc)] 
} 

#----------------------------------------------------------------------- 
# *** PseudoInvers of a Matrix *** 
# using singular value decomposition 
# 
pinv <- function (A) 
{ 
    s <- svd(A) 
    # D <- diag(s$d); Dinv <- diag(1/s$d) 
    # U <- s$u; V <- s$v 
    # A = U D V' 
    # X = V Dinv U' 
    s$v %*% diag(1/s$d) %*% t(s$u) 
} 
#----------------------------------------------------------------------- 


pinv <- function (A) 
{ 
    s <- svd(A) 
     D <- diag(s$d); Dinv <- diag(1/s$d) 
     U <- s$u; V <- s$v 
     A = U D V' 
    X = V Dinv U' 
    s$v %*% diag(1/s$d) %*% t(s$u) 
} 

savgol <- function(T, fl, forder, dorder) { 
    m<-length(T) 
    dorder<-(dorder + 1)

    fc <- (fl-1)/2 
    X <- outer(-fc:fc, 0:forder, FUN="^")
    Y <- pinv(X); 

    T2 <- convolve(T, rev(Y[dorder,]), type="o") 
    T2 <- T2[(fc+1):(length(T2)-fc)] 
}