# Clear the Workspace
	rm(list=ls())

############################################################################
# Define function for smoothing to the maxima
	smoothmax <- function(y, num) {
		# define vector
		val <- c(1:length(y))*0
		# use max within a subset of y for each position in val
		# notice that the last num values of 'val' are left at zero
		for(i in 1:(length(val)-num)) {
			val[i] <- max(y[i:(i+num)])
		}
		val
	}

############################################################################
# Define a function for finding the location where a threshold is passed
	threshold <- function(y, th, start, dir, above) {
		# if dir is TRUE, then go left to right (FALSE is right to left)
		len <- length(y)
		i <- start

		if(dir) {
			# increment right
			if(above) {
				while(y[i] < th) {
					i <- i + 1
					if(i >= len) {
						break 
					}
				}
			} else {
				while(y[i] > th) {
					i <- i + 1
					if(i >= len) {
						break 
					}
				}
			}
		} else {
			# increment left
			if(above) {
				while(y[i] < th) {
					i <- i - 1
					if(i <= 1) {
						break
					}
				}
			} else {
				while(y[i] > th) {
					i <- i - 1
					if(i <= 1) {
						break
					}
				}
			}
		}
		i
	}

#############################################################################
# Load in the audio file
	audio = scan( file="Pulse2k.txt")

# Set up criteria for smoothing and analyzing the data
	mx <- max( audio )
	len <- length( audio )
	smooth <- 12		# At 44.1 kHz, rectified 2 kHz signal repeats every 11 points
	rate <- 44.1
	numstd <- 4			# Number of stdev to determine threshold

# Rectify the signal, smooth it and downsample so each point is 1 millisecond 
	smoothdata <- smoothmax( abs(audio), smooth)
	reducedata <- smoothdata[rate*1:(len/rate-1)]
	plot( reducedata )
	savePlot(filename="Envelope2k.pdf", type="pdf")

#############################################################################
# Find the low end of the growth curve
	start <- 1
	confint <- numstd*sd(reducedata[start:(start+100)])
	if(confint <= 0) {
		confint <- 1
	}
	thresh <- mean(reducedata[start:(start+100)]) + confint
	t1g <- threshold(reducedata, thresh, start, TRUE, TRUE)

# Find the high end of the growth curve
	start <- length(reducedata)/2
	confint <- numstd*sd(reducedata[(start+100):start])
	if(confint <= 0) {
		confint <- 1
	}
	thresh <- mean(reducedata[(start+100):start]) - confint
	t2g <- threshold(reducedata, thresh, start, FALSE, FALSE)

############################################################################
# Find the low end of the decay curve
	start <- length(reducedata)/2
	confint <- numstd*sd(reducedata[(start+100):start])
	if(confint <= 0) {
		confint <- 1
	}
	thresh <- mean(reducedata[(start+100):start]) - confint
	t1d <- threshold(reducedata, thresh, start, TRUE, FALSE)

# Find the high end of the decay curve
	start <- length(reducedata)
	confint <- numstd*sd(reducedata[(start-100):start])
	if(confint <= 0) {
		confint <- 1
	}
	thresh <- mean(reducedata[(start-100):start]) + confint
	t2d <- threshold(reducedata, thresh, start, FALSE, TRUE)

############################################################################
# Plot the growth curve
	growth <- reducedata[t1g:t2g]
	plot(growth, xlab="Time (ms)", ylab="Amplitude", pch="*")
	title("Growth Curve for 2 kHz Signal")
	savePlot(filename="GrowthCurve2k.pdf", type="pdf")

# Determine the growth rate time constant
	transgrowth <- log(1-growth/(mx+1))
	timegrowth <- 1:length(transgrowth)
	plot(timegrowth ~ transgrowth, xlab="Ln(Amplitude)", ylab="Time (ms)", pch="*")
	growthfit <- lm(timegrowth ~ transgrowth)
	abline(growthfit, col=2, lwd=2)
	tau <- growthfit$coefficients[2]
	text(mean(transgrowth),max(timegrowth),paste("tau = ", round(tau, digits=2), " ms"))
	title("Semi-log Growth Curve 2 kHz/Time Constant")
	savePlot(filename="GrowthFit2k.pdf", type="pdf")

###########################################################################
# Plot the decay curve
	decay <- reducedata[t1d:t2d]
	plot(decay, xlab="Time (ms)", ylab="Amplitude", pch="*")
	title("Decay Curve for 2 kHz Signal")
	savePlot(filename="DecayCurve2k.pdf", type="pdf")

# Determine the decay rate time constant
	transdecay <- log(1-decay/(mx+1))
	timedecay <- 1:length(transdecay)
	plot(timedecay ~ transdecay, xlab="Ln(Amplitude)", ylab="Time (ms)", pch="*")
	decayfit <- lm(timedecay ~ transdecay)
	abline(decayfit, col=2, lwd=2)
	tau <- decayfit$coefficients[2]
	text(mean(transdecay),max(timedecay),paste("tau = ",round(tau, digits=2)," ms"))
	title("Semi-log Decay Curve 2 kHz/Time Constant")
	savePlot(filename="DecayFit2k.pdf", type="pdf")

############################################################################
# Refine the analysis of the decay curve
# Clip by hand the window size
	low <- 7000
	high <- 9000
# Plot the decay curve
	decay <- reducedata[low:high]
	plot(decay, xlab="Time (ms)", ylab="Amplitude", pch="*")
	title("Decay Curve for 2 kHz Signal")
	savePlot(filename="DecayCurve2k.pdf", type="pdf")

# Determine the decay rate time constant
#	transdecay <- log(1-decay/(mx+1))
	transdecay <- log(decay/(mx+1))
	timedecay <- 1:length(transdecay)
	plot(timedecay ~ transdecay, xlab="Ln(Amplitude)", ylab="Time (ms)", pch="*")
	decayfit <- lm(timedecay ~ transdecay)
	abline(decayfit, col=2, lwd=2)
	tau <- decayfit$coefficients[2]
	text(mean(transdecay),max(timedecay),paste("tau = ",round(tau, digits=2)," ms"))
	title("Semi-log Decay Curve 2 kHz/Time Constant")
	savePlot(filename="DecayFit2k.pdf", type="pdf")