# Reading of SPE files, produced by Princeton Instruments spectrometers
# File format version 2.5 (Sept. 2002)
# C. Beleites
# R. Kiselev
# July 2015
##' Import WinSpec SPE file
##'
##' Import function for WinSpec SPE files (file version 2.5). The calibration
##' data (polynome and calibration data pairs) for x-axis are automatically
##' read and applied to the spectra. Note that the y-calibration data structure
##' is not extracted from the file since it is not saved there by WinSpec and is
##' always empty.
##'
##' @param filename Name of the SPE file to read data from
##' @param xaxis Units of x-axis, e.g. \emph{"file"}, \emph{"px"},
##' \emph{"nm"}, \emph{"energy"}, \emph{"raman"}, \emph{...}
##' \code{read.spe} function automatically checks if the x-calibration data are
##' available and uses them (if possible) to reconstruct the xaxis
##' in the selected units.
##' @param acc2avg whether to divide the actual data set by the number of
##' accumulations, thus transforming \emph{accumulated} spectra to
##' \emph{averaged} spectra. WinSpec does not do this automatically, so the
##' spectral intensity is always proportional to the number of accumulations.
##' The flag \code{@@data$averaged} is automatically set to \code{TRUE}.
##' @param cts_sec whether to divide the actual data set by the exposure time,
##' thus going to count per second unit.
##' @param keys.hdr2data Which metadata from the file header should be saved to
##' the \code{Data} slot of a newly created hyperSpec object
##'
##' @return hyperSpec object
##'
##' @rdname read.spe
##'
##' @author R. Kiselev, C. Beleites
##' @export
read.spe <- function(filename, xaxis="file", acc2avg=F, cts_sec=F,
keys.hdr2data=c("exposure_sec",
"LaserWavelen",
"accumulCount",
"numFrames",
"darkSubtracted")){
hdr <- read.spe.header(filename)
# Read the part of file that contains actual experimental data
raw.data <- readBin(filename, "raw", file.info(filename)$size, 1)[- (1:4100)]
# Convert raw spectral data according to the datatype defined in the header
spc <- switch(hdr$datatype + 1,
readBin(raw.data, "double", length(raw.data)/4, 4), # float
readBin(raw.data, "integer", length(raw.data)/4, 4, signed=TRUE), # long
readBin(raw.data, "integer", length(raw.data)/2, 2, signed=TRUE), # int
readBin(raw.data, "integer", length(raw.data)/2, 2, signed=FALSE) # uint
)
# Create a structured data.frame that will accomodate spectral data
dim(spc) <- c(hdr$xdim, hdr$ydim * hdr$numFrames)
extra_data <- data.frame (
px.y = rep(seq_len(hdr$ydim), hdr$numFrames),
frame = rep(seq_len(hdr$numFrames), each=hdr$ydim)
)
# Extract selected items from the header. They will go to a new hyperSpec object
hdr2data <- hdr[keys.hdr2data]
if (length (hdr2data > 0))
extra_data <- cbind (extra_data, hdr2data)
# Create hyperSpec object
spc <- new("hyperSpec", spc=t(spc), data=extra_data)
# Modify hyperSpec object
spc@label$spc <- expression("counts")
spc@label$.wavelength <- expression("pixel number")
# Check if we should use display units specified in the SPE file
if (xaxis == "file")
xaxis = .fixunitname(hdr$xCalDisplayUnit)
# Create a new x-axis, if required
xaxis <- .fixunitname(xaxis)
if (xaxis == "px")
return(spc)
if (! hdr$xCalValid)
warning("The calibration is NOT valid")
# Recreate calibration function
polyorder <- hdr$xCalPolyOrder
coeffs <- hdr$xCalPolCoeffs[seq(polyorder + 1)]
vM <- vanderMonde(spc@wavelength, polyorder)
# Check if we have laser wavelength
if (hdr$LaserWavelen < 10)
hdr$LaserWavelen <- NULL
# Perform convertion
spc@wavelength <- wlconv(src=.fixunitname(hdr$xCalPolyUnit),
dst=xaxis,
points=as.numeric(vM %*% coeffs),
laser=hdr$LaserWavelen)
spc@label$.wavelength = switch(xaxis,
nm=expression("Wavelength, nm"),
invcm=expression(tilde(nu) / cm^-1),
ev=expression("Energy / eV"),
freq=expression(nu / THz),
raman=expression(Raman~shift / cm^-1))
if (acc2avg){
spc <- spc / hdr$accumulCount
spc@data$averaged <- T
}
if (cts_sec){
spc <- spc / hdr$exposure_sec
spc@label$spc <- expression("counts / s")
}
## consistent file import behaviour across import functions
.fileio.optional (spc, filename)
}
##' @describeIn read.spe Read only header of a WinSpec SPE file (version 2.5)
##' @return hdr list with \code{key=value} pairs
##' @export
read.spe.header <- function(filename){
# Read the 4100-byte long binary header from the SPE file and parse it
# Load the header
raw.data <- readBin(filename, "raw", 4100, 1)
# Extract some items from the 4100 bytes-long file header
hdr <- list (
hwVersion = readBin(raw.data[1 :2 ], "integer", 1, 2, signed=TRUE ), # uint16
xDimDet = readBin(raw.data[7 :8 ], "integer", 1, 2, signed=FALSE), # uint16
mode = readBin(raw.data[9 :10 ], "integer", 1, 2, signed=TRUE ), # uint16
exposure_sec = readBin(raw.data[11 :14 ], "double", 1, 4), # float32
vChipXDim = readBin(raw.data[15 :16 ], "integer", 1, 2, signed=TRUE ), # int8
vChipYDim = readBin(raw.data[17 :18 ], "integer", 1, 2, signed=TRUE ), # int8
yDimDet = readBin(raw.data[19 :20 ], "integer", 1, 2, signed=FALSE), # uint16
date = readBin(raw.data[21 :30 ], "character", 1, 10 ), # char
detTemperature = readBin(raw.data[37 :40 ], "double", 1, 4), # float32
xdim = readBin(raw.data[43 :44 ], "integer", 1, 2, signed=FALSE), # uint16
shutterMode = readBin(raw.data[51 :52 ], "integer", 1, 2, signed=FALSE), # uint16
specCenterWlNm = readBin(raw.data[73 :76 ], "double", 1, 4), # float32
datatype = readBin(raw.data[109 :110 ], "integer", 1, 2, signed=TRUE ), # int8
darkSubtracted = readBin(raw.data[151 :152 ], "integer", 1, 2, signed=FALSE), # int8
timeLocal = readBin(raw.data[173 :179 ], "character", 1, 7 ), # char
timeUTC = readBin(raw.data[180 :186 ], "character", 1, 7 ), # char
gain = readBin(raw.data[199 :200 ], "integer", 1, 2, signed=FALSE), # uint16
comments = readBin(raw.data[201 :600 ], "character", 1, 400 ), # char
ydim = readBin(raw.data[657 :658 ], "integer", 1, 2, signed=FALSE), # uint16
accumulCount = readBin(raw.data[669 :672 ], "integer", 1, 4), # uint32
readoutTime = readBin(raw.data[673 :676 ], "double", 1, 4), # float32
swVersion = readBin(raw.data[688 :704 ], "character", 1, 16 ), # char
kinTrigMode = readBin(raw.data[725 :726 ], "integer", 1, 2, signed=TRUE ), # int16
expRepeatCount = readBin(raw.data[1419:1422], "integer", 1, 4, signed=TRUE ), # int32
expAccumCount = readBin(raw.data[1423:1426], "integer", 1, 4, signed=TRUE ), # int32
hwAccumFlag = readBin(raw.data[1433:1434], "integer", 1, 2, signed=TRUE ), # int16
cosmicApplied = readBin(raw.data[1439:1440], "integer", 1, 2, signed=TRUE ), # int16
cosmicType = readBin(raw.data[1441:1442], "integer", 1, 2, signed=TRUE ), # int16
numFrames = readBin(raw.data[1447:1450], "integer", 1, 4), # int32
shutterType = readBin(raw.data[1475:1476], "integer", 1, 2, signed=TRUE ), # int16
readoutMode = readBin(raw.data[1481:1482], "integer", 1, 2, signed=TRUE ), # int16
kinWindowSize = readBin(raw.data[1483:1484], "integer", 1, 2, signed=TRUE ), # int16
clkSpeed = readBin(raw.data[1485:1486], "integer", 1, 2, signed=TRUE ), # int16
computerIface = readBin(raw.data[1487:1488], "integer", 1, 2, signed=TRUE ), # int16
# X Calibration Structure
xCalOffset = readBin(raw.data[3001:3008], "double", 1, 8, signed=TRUE ), # float64
xCalFactor = readBin(raw.data[3009:3016], "double", 1, 8, signed=TRUE ), # float64
xCalDisplayUnit= readBin(raw.data[3017 ], "integer", 1, 1, signed=FALSE), # uint8
xCalValid = readBin(raw.data[3099 ], "integer", 1, 1, signed=FALSE), # uint8
xCalInputUnit = readBin(raw.data[3100 ], "integer", 1, 1, signed=FALSE), # uint8
xCalPolyUnit = readBin(raw.data[3101 ], "integer", 1, 1, signed=FALSE), # uint8
xCalPolyOrder = readBin(raw.data[3102 ], "integer", 1, 1, signed=FALSE), # uint8
xCalPointCount = readBin(raw.data[3103 ], "integer", 1, 1, signed=FALSE), # uint8
xCalPxPos = readBin(raw.data[3104:3183], "double", 10, 8, signed=TRUE ), # float64
xCalValues = readBin(raw.data[3184:3263], "double", 10, 8, signed=TRUE ), # float64
xCalPolCoeffs = readBin(raw.data[3264:3311], "double", 6, 8, signed=TRUE ), # float64
LaserWavelen = readBin(raw.data[3312:3319], "double", 1, 8, signed=TRUE ) # float64
)
# Convert magic numbers into human-readable unit strings
spe_units <- c("pixel", "data", "user units", "nm", "cm-1", "Raman shift")
hdr$xCalDisplayUnit <- spe_units[hdr$xCalDisplayUnit]
hdr$xCalInputUnit <- spe_units[hdr$xCalInputUnit]
hdr$xCalPolyUnit <- spe_units[hdr$xCalPolyUnit]
return(hdr)
}
##' @describeIn read.spe Plot the WinSpec SPE file (version 2.5) and show the
##' calibration points stored inside of it (x-axis calibration)
##' @export
spe.showcalpoints <- function(filename, xaxis="file", acc2avg=F, cts_sec=F){
hdr <- read.spe.header(filename)
xaxis <- .fixunitname(xaxis)
# Check if we should use display units specified in the SPE file
if (xaxis == "file")
xaxis <- .fixunitname(hdr$xCalDisplayUnit)
if (xaxis == "px"){
xaxis <- hdr$xCalPolyUnit
warning("Cannot show calibration data in pixels")
}
# Open file, make plot and mark position of all peaks stored inside the file
# in the x-calibration structure
spc <- read.spe(filename, xaxis, acc2avg, cts_sec)
rng <- max(spc) - min(spc)
ylims <- c(min(spc), max(spc) + 0.3*rng)
if (dim(spc@data$spc)[1] > 1)
plot(spc, plot.args=list(ylim=(ylims)), "spcprctl5")
else
plot(spc, plot.args=list(ylim=(ylims)))
title(basename(filename))
if (hdr$xCalPointCount == 0){
warning("No calibration data! Nothing to show")
return("")
}
markpeak(spc, wlconv(src=hdr$xCalInputUnit,
dst=.fixunitname(xaxis),
points=hdr$xCalValues,
laser=hdr$LaserWavelen))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.