# R program for comparing mass spectra

msbin <- function(x, y, bin.width=0.25, x.max=max(x), x.min=min(x)){
# ---------------------------------------------------------------
# Simple binning of an MS/MS spectrum
#
# x             input array of m/z values
# y             input array of intensity values
# bin.width     m/z bin width of output
#
# Value:
# x     =       evenly sampled m/z values for the binned spectrum
# y     =       binned intensity values
# ---------------------------------------------------------------

  x.n <- length(x)

  n.bins <- round( (x.max - x.min) / bin.width ) + 1

  x.bin <- seq(from=x.min, to=x.max, length=n.bins)
  y.bin <- rep(0, n.bins)

  for(i in 1:x.n){
    idx <- round( (x[i] - x.min) / bin.width ) + 1
    y.bin[idx] <- y.bin[idx] + y[i]
  }

  out <- list(x = x.bin, y = y.bin)
  out
}

fftconvolve <- function(y1, y2){
# ----------------------------------------------------------
# Function to perform FFT based convolution of two sequences
# ----------------------------------------------------------

  y1.n <- length(y1)
  y2.n <- length(y2)
 
  n.out <- y1.n + y2.n - 1

  nf <- log(n.out) / log(2.0)
  if (nf == floor(nf)){
        n <- floor(nf)
  }
  else{
        n <- floor(nf) + 1
  }

  n.round <- 2^n

  if (n.round > y1.n){   
    y1 <- c(y1, rep(0.0, n.round - y1.n))
  }
  if (n.round > y2.n){
        y2 <- c(y2, rep(0.0, n.round - y2.n))
  }

  y1.fft <- fft(y1)
  y2.fft <- fft(y2)
     
  y.real <- (Re(y1.fft)*Re(y2.fft) - Im(y1.fft)*Im(y2.fft)) / n.round
  y.imag <- (Im(y1.fft)*Re(y2.fft) + Re(y1.fft)*Im(y2.fft)) / n.round
  y.ifft <- complex(length=n.round, real=y.real, imaginary=y.imag)
 
  y.out <- fft(y.ifft, inverse=T)
 
  out <- Re(y.out[1:n.out])
}

gmg <- function(x, a0, a1, a2, a3){
# -----------------------------------------------------
# Produce a Gaussian Modified Gaussian peak given x and
# a0 (peak area)
# a1 (mean)
# a2 (sigma)
# a3 (skewness)
# -----------------------------------------------------
  y <- a0 * exp(-0.5 * (x-a1)^2/(a3^2 + a2^2)) *
                2*pnorm(a3*(x-a1)/(a2*sqrt(a3^2 + a2^2))) /
                (sqrt(2 * 3.14159) * sqrt(a3^2 + a2^2))
  y
}

peak <- gmg(seq(-5,5,0.1),3,1,.2,0)

dmmp <- as.matrix(read.table("dmmp.dat", sep=","))
dmmp.bin <- msbin( dmmp[,1], dmmp[,2], bin.width=0.1, x.min=20, x.max=200 )
dmmp.bin$y <- dmmp.bin$y/max(dmmp.bin$y)
dmmp.m<-fftconvolve(dmmp.bin$y, peak)
dmmp.mz<-data.frame(dmmp.bin$x,dmmp.m[61:(length(dmmp.m)-40)])
names(dmmp.mz)<-c("mz","I")
names(dmmp.bin)<-c("mz","I")
dmmp.mz$I <- dmmp.mz$I/max(dmmp.mz$I)

demp <- as.matrix(read.table("demp.dat", sep=","))
demp.bin <- msbin( demp[,1], demp[,2], bin.width=0.1, x.min=20, x.max=200 )
demp.bin$y <- demp.bin$y/max(demp.bin$y)
demp.m<-fftconvolve(demp.bin$y, peak)
demp.mz<-data.frame(demp.bin$x,demp.m[61:(length(demp.m)-40)])
names(demp.mz)<-c("mz","I")
names(demp.bin)<-c("mz","I")
demp.mz$I <- demp.mz$I/max(demp.mz$I)

dmmp.filt<-dmmp.mz

zeros<-rep(0,length(dmmp.mz$mz))

dmmp.scan<-c(zeros,dmmp.mz[,2],zeros)
demp.scan<-c(zeros,demp.mz[,2],zeros)

par(mfcol=c(4,2))

for( i in c(0.01, 0.25))
{
	dmmp.noise<-dmmp.scan + rnorm(length(dmmp.scan),sd=i)
	demp.noise<-demp.scan + rnorm(length(demp.scan),sd=i)

	xr <- c(1801,1801*2)
	yr <- c(0,max(dmmp.noise))

	plot(dmmp.noise, ylim=yr, xlim=xr, type="l")
	title("DMMP Spectrum")
	dmmp.clean<-fftconvolve(dmmp.noise,rev(dmmp.filt$I))
	fyr <- c(0,max(dmmp.clean))
	plot(dmmp.clean, ylim=fyr, xlim=c(3500,3750), type="l")
	title("DMMP from DMMP filter")

	plot(demp.noise,ylim=yr,  xlim=xr, type="l")
	title("DEMP Spectrum")
	demp.clean<-fftconvolve(demp.noise,rev(dmmp.filt$I))
	plot(demp.clean, ylim=fyr, xlim=c(3500,3750), type="l")
	title("DEMP from DMMP filter")
}