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.

# FFT and derivatives

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')      

# Noisy signals

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)

## Filtering

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.

## Noise cancellation by Averaging

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

# Image Processing

Refs:

Some useful functions:

# to install: source("http://bioconductor.org/biocLite.R"); biocLite("EBImage")
library("EBImage")

# plot a picture
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:
show_pic(clown)