fourier01_03.r

#select peaks in spectrum to filter out
#click on each peak, then end the locator
peaks <- data.frame(locator())

#function to create a filter for circle shaped peaks in the spectrum
#dimensions: dimensions of the spectrum matrix
#x: x position of the centre of the circle to be filtered out
#y: y position of ...
#r: radius of circle to be filtered out of the spectrum
#gaussian: flag that indicates whether the circle should have a Gaussian density or not
#          when set to TRUE, smoother results are obtained.
filter_circle <- function(dimensions, x, y, r = 2, gaussian = T)
{
  result   <- matrix(1, dimensions[1], dimensions[2])
  distance <- sqrt((row(result) - x)^2 + (col(result) - y)^2)
  if (gaussian)
  {
    result <- dnorm(distance, sd = r*2)
    result <- 1 - result/max(result)
  } else
  {
    result[distance < r] <- 0
  }
  return(result)
}

jpeg("fourier_spectrum.jpg", nrow(magni), ncol(magni))
par(mar = c(0,0,0,0), oma = c(0,0,0,0))
image(1:nrow(magni), 1:ncol(magni), log(magni + 1), col = grey((0:255)/255), bty = "n", xaxt = "n", yaxt = "n")
points(y~x, peaks, col = "red", cex = 3)
dev.off()

#create filter
filt <- by(peaks, 1:nrow(peaks), function(x) filter_circle(dim(magni), x$x, x$y, 14))

#function that returns the product of all elements in a list
list_prod <- function(x)
{
  if (length(x) > 1)
  {
    result      <- x[2:length(x)]
    result[[1]] <- result[[1]]*x[[1]]
    return(list_prod(result))
  } else return(x)
}
filt <- list_prod(filt)[[1]]

#filter out the selected peaks in the spectrum:
magni_filter <- magni*filt