# R code for Fourier filtering

# load SignalNoise.RData and examine the four objects within:
#   noise: random noise from population with a mean of 0 and a 
#          standard deviation of 10
#   time: an index representing increments of time
#   signal: the pure signal, which consists of three Gaussian 
#           peaks at times of 256, 513, and 641
#   sig.tot: the sum of signal and time, or the noisy signal
load("SignalNoise.RData")
mean(noise)
sd(noise)

# first, look at the help file; note that fft() completes the
# Fourier transform from the original domain (time or frequency)
# to the other domain (frequency or time); the option inverse =
# TRUE reverses the process; note that when z is a vector (which
# is true for us), the result of fft() in unnormalized
help(fft)

# next, let's create a simple test vector to see how this works
test = c(1, 2, 3, 4)
test

# complete the FT to move from one domain to the other; note that
# the result of the FT is a complex number
ft.test = fft(test)
ft.test

# complete the inverse FT to return to the original domain; note
# that each individual value is now four times larger; this is
# because the result is returned without normalization
ift.test = fft(ft.test, inverse = TRUE)
ift.test

# let's fix the normalization problem by dividing by the length
# of our vector
ift.test = fft(ft.test, inverse = TRUE)/length(test)
ift.test

# we can extract the real portion of the result using Re(), which
# returns the real portion of a complex number
Re(ift.test)

# apply fft to our noisy signal
sim.ft = fft(sig.tot)
head(sim.ft)

# plot the real components of the fourier transformed data
plot(time, Re(sim.ft), type = "l", col = "blue")

# note the symmetry of the resulting plot; the first half of the
# data (512 points) contains the same information as the last
# half of the data; let's look more closely at the first half
plot(time[1:512], Re(sim.ft)[1:512], type = "l", col = "blue")

# the points at the beginning are dominated by the signal, which
# is why there is a systematic variation in height with time; the
# remaining points are dominated by noise, which is why the
# variation in height with time is random; to filter out the
# noise, we set the values of these points to zero; if we do this
# for poionts 101-512, then we also need to do it for points
# 513-924, leaving alone 100 points at either end to give us the
# signal
sim.ft[101:924] = 0 + 0i

# here is the fourier transformed data after zeroing out the noise
plot(time[1:512], Re(sim.ft)[1:512], type = "l", col = "blue")
plot(time, Re(sim.ft), type = "l", col = "blue")

# now we take the inverse fourier transform
sim.ift = fft(sim.ft, inverse = TRUE)/length(sim.ft)

# show plots to compare original data and filtered data
old.par = par(mfrow = c(2, 1))
plot(time, sig.tot, type = "l", col = "blue")
plot(time, Re(sim.ift), type = "l", col = "blue")
par(old.par)

# finally, the signal-to-noise ratio before and after filtering
sd.before = sd(sig.tot[800:1000])
m.before = mean(sig.tot[254:258])
sn.before = m.before/sd.before

sd.after = sd(Re(sim.ift[800:1000]))
m.after = mean(Re(sim.ift[254:258]))
sn.after = m.after/sd.after

ff_df = data.frame(c("before", "after"), 
                   c(m.before, m.after), 
                   c(sd.before, sd.after), 
                   c(sn.before, sn.after))

colnames(ff_df) = c("","mean signal", 
                    "std dev noise", "S/N")

ff_df