#  Fishing for peaks

smooth <- function( y )
   {
   	n <- length( y )
   	( y[1:(n-2)] + y[2:(n-1)] + y[3:n] ) / 3
   }
   
dat <- read.table("T08030205.txt",skip=1)

sig <- dat$V9                 # get channel 1 signal
noise <- dat$V11              # get channel 1 noise
plot( sig, type="l", xlim=c( 9000, 10000 ) )    # take a look at the traces  
plot( noise, type="l", xlim=c( 9000, 10000 ) )  #
plot( noise, sig, pch=20 )                      # are they correlated?

sp <- smooth( sig )           # smooth the data a bit
np <- smooth( noise )

plot( sp, type="l", xlim=c(9000,10000) )   
plot( np, type="l", xlim=c(9000,10000) )
plot( sp - np, type="l", xlim=c(9000,10000) )

yp <- sp - np
yp <- yp - min( yp )
plot( sp - np, type="l", xlim=c(9000,10000) )

m <- (yp > 8) * 1             #  1 where above threshold of 8, 0 where not
n <- length(m)

dm <- m[2:n] - m[1:(n-1)]     #  incremental differences betweed adjacent
                              #  ones and zeros1 for start, -1 at stop

m[1:100]
dm[1:100]
      
starts <- (1:n)[ dm == 1 ]    # (1:n) is the vector 1 2 3 4 ... n so
stops <- (1:n)[ dm == -1 ]    # indexing with boolean vector yields positions

length(starts)                # check that lengths are the same
length(stops)

hist(stops-starts)                 # histogram of lengths above threshold
hist(stops-starts,breaks=0:450 )
sort(stops-starts)

# pt is a matrix of starts, stops and lengths above threshold

pt <- cbind( starts, stops, len = (stops-starts) + 1 )


getOption("digits")           # what is output precision?
options(digits = 3)           # change precision

pt

sort(pt[,3])                  # sort column 3 (lengths) 
order(pt[,3])                 # get indices of sorted 3
pt[order(pt[,3]),]            # order entire pt TABLE by column 3

plot(yp,type="l",xlim=c(25000,27000))   # zoom in on regions
plot(yp,type="l",xlim=c(29000,30000))

#  apply to rows of matrix pt the given function
#          1 means rows, 2 means columns
#
#   this expression returns, for each row in pt, the yp vales for all positions 
#   between the start and stop
apply( pt, 1, function( row ) yp[ (row[1]+1) : row[2] ] )
#
#   this returns, for each row, the total integrated signal of yp between
#   each start and stop (area under the curve above threshold)
apply( pt, 1, function( row ) sum( yp[ (row[1]+1) : row[2] ] ) )

# save result as vector
areas <- apply( pt, 1, function( row ) sum( yp[ (row[1]+1) : row[2] ] ) )
#
# add the area vector as a four column to pt
#
pt <- cbind( pt, areas )

# 
# zoom in, look at effects of more smoothing
#
plot( yp, type="l", xlim=c(29000,30000) )
ypp <- smooth( smooth( yp ) )
plot( ypp, type="l", xlim=c(29000,30000) )
points( yp, type="l", col="green" )
points( ypp, type="l" )
ypppp <- smooth( smooth( smooth( smooth( yp ) ) ) )
points( ypppp, type="l", col="red" )