# ---------------------------------------------------------------------
# Program: FFT2.S  
#  Author: Steven M. Boker
#    Date: Tue Oct 14 13:42:50 EST 2003
#
# This simulates two composite and plots an FFT of a file from disk.
# 
#
# ---------------------------------------------------------------------
# Revision History
#   Steven M. Boker -- Tue Oct 14 13:42:59 EST 2003
#      Created FFT2.S
#
# ---------------------------------------------------------------------


# ---------------------------------------------------------------------
# Read in a data file

theDataList <- scan("sine1.dat", list(x=0), skip=0)

maxSamples <- length(theDataList$x)

maxSamples

summary(theDataList$x)

# ---------------------------------------------------------------------
# Plot the data

thePlotData <- theDataList$x
tepsfile <- paste("Sine1TimeSeries.eps",sep="")
tmain <- paste("Sine1 Time Series", sep="")
postscript(tepsfile,height=6.4,width=7.5,horizontal=T)
plot(c(0, maxSamples), c(-2,2), xlab="Samples", ylab="Displacement", 
    main=tmain, type="n")
lines(seq(1, maxSamples),thePlotData[1:maxSamples], type="l")
dev.off()

# ---------------------------------------------------------------------
# Calculate a periodogram using the built in fft function

# -----------------------------
# Define the periodogram function

pgram <-function(z) {
	n <- length(z)
	(Mod(fft(z)^2/(2*pi*n)))[1:(n %/% 2 + 1)]
}

# -----------------------------
# Calculate the periodogram of tSine1

tSine1Pgram <- pgram(theDataList$x)

# -----------------------------
# Plot the full Frequency Power Spectrum for theDataList$x

thePlotData <- (tSine1Pgram / sum(tSine1Pgram)) * 100
theFreq <- rep(0,maxSamples/2)
for(j in 1:(maxSamples/2)) {
	theFreq[j] <- (j/maxSamples) * (2*3.14159)
}
tepsfile <- paste("Sine1Pgram.eps",sep="")
tmain <- paste("Frequency Power Spectrum of Sine1", sep="")
postscript(tepsfile,height=6.4,width=7.5,horizontal=T)
plot(c(0, theFreq[maxSamples/8]), c(0,50), xlab="Frequency", 
      ylab="Percent of Total Power", main=tmain, type="n")
lines(theFreq,thePlotData[1:(maxSamples/2)], type="l")
dev.off()


# ---------------------------------------------------------------------
# End the program

q()