# ---------------------------------------------------------------------
# Program: FFT1.S  
#  Author: Steven M. Boker
#    Date: Thu Mar  4 10:12:51 EST 1999
#
# This simulates composites and then calculates and plots an FFT.
# 
#
# ---------------------------------------------------------------------
# Revision History
#   Steven M. Boker -- Thu Mar  4 10:12:54 EST 1999
#      Created FFT1.S
#
# ---------------------------------------------------------------------


# ---------------------------------------------------------------------
# Simulate two composite sine waves

maxSamples <- 512

# -----------------------------
# Select two frequencies and phase offset

theFreq1 <- .5
theFreq2 <- .1
phaseOffset <- 20

# -----------------------------
# The first data vector is a composite of the two frequencies

tSine1 <- sin(theFreq1*(1:(maxSamples+phaseOffset))) + 
          sin(theFreq2*(1:(maxSamples+phaseOffset)))

# -----------------------------
# The second data vector is dependent on the first time series


tSine2 <- 0.1 * tSine1[(1+phaseOffset):(maxSamples+phaseOffset)] + 
          rnorm((maxSamples), mean=0, sd=.2)

tSine1 <- tSine1[1:maxSamples]

# ---------------------------------------------------------------------
# Plot the two composite sine waves

thePlotData <- tSine1
tepsfile <- paste("Sine1TimeSeries.eps",sep="")
tmain <- paste("Composite Sine Wave 1", 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()

thePlotData <- tSine2
tepsfile <- paste("Sine2TimeSeries.eps",sep="")
tmain <- paste("Composite Sine Wave 2", 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(tSine1)

# -----------------------------
# Plot the full Frequency Power Spectrum for tSine1

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("Sine1PgramFull.eps",sep="")
tmain <- paste("Frequency Power Spectrum of Sine Wave 1", sep="")
postscript(tepsfile,height=6.4,width=7.5,horizontal=T)
plot(c(0, theFreq[maxSamples/2]), c(0,50), xlab="Frequency", 
      ylab="Percent of Total Power", main=tmain, type="n")
lines(theFreq,thePlotData[1:(maxSamples/2)], type="l")
dev.off()


# -----------------------------
# Calculate the periodogram of tSine2

tSine2Pgram <- pgram(tSine2)

# -----------------------------
# Plot the full Frequency Power Spectrum for tSine2

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

# -----------------------------
# Calculate the coherence between sine wave 1 and sine wave 2

cor(tSine1Pgram, tSine2Pgram)



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

q()