# 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