Refs:
These page presents some egs taken mostly from the mentioned MOOC.
L <- 20; # interval of interest
n <- 128; # sampling rate
x = head(seq(-L/2,L/2,len=n+1), -1); # the sampling rate was taken on regular intervals
u = exp(-x ^ 2); # some signal (in this case, a gaussian)
plot(x,u,type='l', xlab="samples", ylab="signal")
ut <- fft(u) # compute the signal's fourier transform
k <- (2*pi/L) * c(seq(0,n/2-1),seq(-n/2,-1)); # rescaling frequencies from [0,2pi] to [-L/2,L/2]
library(pracma) # sech, tanh, fftshift, ifft
par(mfrow=c(2,2))
plot(x,ut,type='l',main="signal spectra", xaxt='n',xlab=NA, ylab="strength")
# fft rearranges the outputs of fft by moving the zero-frequency component
# to the center of the array
plot(fftshift(k), fftshift(ut), type='l', main="spectra centered at frequency zero", xlab="frequencies", ylab="strength")
# take the norm, here the phase doesn't matter
plot(fftshift(k), abs(fftshift(ut)),type='l', main="norm of spectra", xlab="frequencies", ylab="strength")
Yeap, the FFT of a Gaussian is a Gaussian.
The FFt can be used to approximate derivatives.
Given function \(f\) and its n-th derivative \(f^{(n)}\), \[\text{fft}( f^{(n)} ) \approx (ik)^n \text{fft}( f )\] where \(i\) is \(\sqrt{-1}\) and \(k\) the vector of frequencies
u = sech(x) # a function
ud = -sech(x) * tanh(x) # its first derivative
u2d = sech(x) - 2*sech(x)^3 # and its 2nd derivative
ut = fft(u)
uds = ifft( (1i*k) *ut) # ifft: Inverse fast Fourier transform
u2ds = ifft( (1i*k)^2*ut)
# draw first derivative and the fourier approximation
plot(x,ud,type='l',col='black', xlab="x", ylab="df/dx")
points(x,uds,col='red')
# draw second derivative and the fourier approximation
plot(x,u2d,type='l',col='black', xlab="x", ylab="d^2f/dx^2")
points(x,u2ds,col='red')
This is the original signal:
T <- 30; # signal lasted 30 seconds and was sampled 512 times
n <- 512;
t <- head(seq(-T/2,T/2,len=n+1),-1); # the sampling rate was taken on regular intervals
k <- (2*pi/T) * c(seq(0,n/2-1),seq(-n/2,-1)) # rescaling for interval T
ks <- fftshift(k) # center frequencies around zero
u <- sech(t); # some signal (in this case, a hiperbolic cos)
ut <- fft(u)
par(mfrow=c(2,1))
plot(t,u,type='l',main="signal in time", xlab="time samples", ylab="signal")
plot(ks,abs(fftshift(ut))/max(abs(fftshift(ut))),type='l',main="signal in frequency", xlab="frequencies", ylab="strength")
Now, let’s add noise to the frequency spectra:
noise <- 20
utn <- ut + noise*(rnorm(n) + 1i*rnorm(n))
It’s important to add imaginary noise, which means noise has also a phase.
Adding only real noise will make it symmetric in the time domain; adding only imaginary noise will make it anti-symmetric.
un = ifft(utn); # the noisy signal from the noisy spectra
par(mfrow=c(2,1))
plot(ks,abs(fftshift(ut))/max(abs(fftshift(ut))),type='l')
points(ks,abs(fftshift(utn))/max(abs(fftshift(utn))),col='magenta',type='l')
legend("topleft",c("real spectra", "noisy spectra"), col=c("black","magenta"),lty=1)
plot(t,u,type='l',ylim=c(0,3)) # real signal
points(t,abs(un),col='magenta',type='l') # noisy signal
legend("topleft",c("real signal", "noisy signal"), col=c("black","magenta"),lty=1)
Let’s assume we just have the noisy data
If we are interesting just in certain frequencies (like in Radar) we can filter the signal (here ‘utn’) to have just the interval we want
Assume we are interested in frequency zero, and are using a gaussian filter:
filter = exp(-k^2) # gaussian filter
utnf = filter * utn # the filterted frequencies
unf = ifft(utnf); # and we rebuild the filtered signal
par(mfrow=c(2,1))
plot(t,abs(unf),col='magenta',type='l')
# we also include a decision line (herein at 0.5) which is the threshold
# in where we decide that it's a relevant signal
abline(h=0.5,col="red",lty=2)
plot(ks,abs(fftshift(utnf))/max(abs(fftshift(utnf))),col='magenta',type='l')
Let’s try to apply the same process at another place, say at -10:
mu = -10
filter = exp(-(k-mu)^2) # gaussian filter
utnf = filter * utn # the filterted frequencies
unf = ifft(utnf); # and we rebuild the filtered signal
par(mfrow=c(2,1))
plot(t,abs(unf),col='magenta',type='l',ylim=c(0,0.6))
abline(h=0.5,col="red",lty=2)
plot(ks,abs(fftshift(utnf))/max(abs(fftshift(utnf))),col='magenta',type='l')
This happens because the signal is coherent, ie, it is in phase, reinforcing itself. The noise on the other hand is out of phase.
If we don’t know the frequency of the signal and we can sample several times, averaging is a powerful way to remove noise, since noise will cancel itself out, while the signal continues. We could average in time/space or in frequency. In frequency is more robust because the signal can be moving and not be in the same place on each sampling, while the object (in principle) keeps signalling at the same frequencies.
avg <- rep(0, n)
realizations <- 30
for (j in 1:realizations) {
avg <- avg + ut + noise*(rnorm(n) + 1i*rnorm(n))
}
avg <- abs(fftshift(avg))/realizations
plot(ks,avg, type='l') # lo and behold
Refs:
Some useful functions:
# to install: source("http://bioconductor.org/biocLite.R"); biocLite("EBImage")
library("EBImage")
# plot a picture
show_pic <- function(pic, adjust=FALSE) {
dims <- dim(pic)
plot(c(0, dims[1]), c(0, dims[2]), type='n', xlab="", ylab="")
rasterImage(pic,0,0,dims[1],dims[2]) # present image
}
# log transformation
# http://homepages.inf.ed.ac.uk/rbf/HIPR2/pixlog.htm
log_transf <- function(pic) {
c <- 1/(log(1+max(abs(pic)))) # pixels in [0,1]
c * log(1+abs(pic))
}
# compute the abs of a matrix and normalizes it to [0,1]
abs_norm <- function(pic) {
nm <- abs(pic)
nm/max(nm)
}
# place the dc-component in the middle (like a 2D fftshift)
# if m is a square 2^n matrix,
# center_matrix(center_matrix(m)) == m
center_matrix <- function(m) {
dims <- dim(m)
m1 <- cbind(m[,(dims[2]/2+1):dims[2]],m[,1:(dims[2]/2)])
rbind(m1[(dims[1]/2+1):dims[1],],m1[1:(dims[1]/2),])
}
# an eg:
clown <- readImage("cln1.png")
show_pic(clown)