#
# calculate isotopic distributions of molecules using the FFT
#
# (c) Magnus Palmblad, 1999
#
# translated from MATLAB to R (by Magnus Palmblad) 2010
#
MAX_ELEMENTS <- 5; MAX_MASS <- 2^13
M <- c(378,234,65,75,6) # empirical formula, e.g. bovine insulin
A<-rep(0,times=MAX_ELEMENTS*MAX_MASS); dim(A)<-c(MAX_ELEMENTS,MAX_MASS)
A[1,2:3] <- c(0.9998443,0.0001557)
A[2,13:14] <- c(0.98889,0.01111)
A[3,15:16] <- c(0.99634,0.00366)
A[4,17:19] <- c(0.997628,0.000372,0.002000)
A[5,33:37] <- c(0.95018,0.00750,0.04215,0,0.00017)
tA <- mvfft(t(A)) # FFT along each element's isotopic distribution
ptA <- t(rep(1,MAX_MASS))
for(i in 1:MAX_ELEMENTS) ptA <- ptA*(tA[,i]^M[i]) # multiply transforms (elementwise)
riptA <- Re(mvfft(t(ptA), inverse = TRUE)) # inverse FFT to get convolutions
plot(0:(MAX_MASS-1), riptA/sum(riptA), type = "h", xlab="integer mass", ylab="frequency") # plot of isotopic distribution
# Addendum: To zoom in around the isotopic distribution, use something like this:
#
# plot_from <- which.max(riptA)-10
# plot_to <- which.max(riptA)+20
# plot(plot_from:plot_to, riptA[(plot_from+1):(plot_to+1)]/sum(riptA[(plot_from+1):(plot_to+1)]), type = "h", xlab="integer mass", ylab="frequency")
#