R/readEDF.R
In gjmvanboxtel/eegr: Electrophysiological data analysis
# readEDF.R - Function to read an EDF/BDF file into an R data frame
# Copyright (C) 2012, 2019 Geert van Boxtel,
# Tilburg University, G.J.M.vBoxtel@tilburguniversity.edu
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
# See also: http://www.gnu.org/licenses/gpl-2.0.txt
#
# Version history:
# 20120207 GvB Original version
# 20190701 GvB complete overhaul for eegr v0.1.0
# 20201104 GvB link to gsignal instead of signal package
# 20230516 GvB bugfixes in calculating record size and up/downsampling
#
#-----------------------------------------------------------------------------------------------
#' readEDF
#'
#' Read data from a European Data Format (EDF, EDF+, BDF, BDF+) file.
#'
#' The function readEDF reads files containing electrophysiological data stored in the European Data
#' Format (EDF) or its variants.
#'
#' EDF was published in 1992 and stores multichannel data, allowing different sample rates for each signal.
#' Internally it includes a header and one or more data records. The header contains some general information
#' (patient identification, start time...) and technical specs of each signal (calibration, sampling rate, filtering, ...),
#' coded as ASCII characters. The data records contain samples as little-endian 16-bit integers.
#'
#' EDF+ was published in 2003 and is largely compatible to EDF, but EDF+ files also allow coding discontinuous recordings
#' as well as annotations, stimuli and events in UTF-8 format.
#'
#' Biosemi developed an EEG acquisition system with 24-bit A/D converters, which could not be stored in the 16-bit
#' EDF format. They introduced a 24-bit variant of EDF called BDF (Biosemi Data Format), and a BDF+ variant is sometimes also
#' found.
#'
#' The function \code{readEDF} can read all of these variants and produce different output that can be used for further processing
#' by function in the \code{eegr} package.
#'
#' @references Kemp B., Värri, A., Rosa, A.C., Nielsen, K.D. and Gade, J. (1992). A simple format for exchange of
#' digitized polygraphic recordings. \emph{Electroencephalogr Clin Neurophysiol}, \emph{82(5)}, 391-393.
#' @references Kemp, B. and Olivan, J. (2003). European data format ’plus’ (EDF+), an EDF alike standard format
#' for the exchange of physiological data. \emph{Clin Neurophysiol}, \emph{114(9)}, 1755-1761.
#' @references European Data Format. \url{https://www.edfplus.info/}
#' @references Which file format does Biosemi use? \url{https://www.biosemi.com/faq/file_format.htm}
#'
#' @author Geert van Boxtel, \email{G.J.M.vanBoxtel@@gmail.com}
#'
#' @param file Character string. The name of the file which the data are to be read from. \code{file}
#' can also be a complete URL. (For the supported URL schemes,see the ‘URLs’ section of the help for url.)
#' @param records Numeric or \code{'all'}. In EDF, signals are stored in 'records', usually containing 0.5 or 1
#' second of data. This function either reads all records (default: \code{all}), or a subset of the records, thus
#' limiting the time interval of the signals. \code{records} can be a range (e.g., \code{1:50}) or a random sequence
#' of records (e.g., \code{c(1,5,12,3)}). In the latter case, the records are returned in the specified order (only
#' useful for discontinuous data). If \code{length(records) == 0} or \code{is.na(records)}, then no records are read
#' and only the file header is returned.
#' @param signals Numeric, character string, or \code{'all'}. Specifies the signals to read from the input file.
#' \code{signals} can be a range (e.g., \code{1:21}) a random sequence of signals (e.g., \code{c(1,5,12,3)}),
#' or a sequence of signal labels (e.g., \code{c("F3", "Fz", "F4")}. In the latter two cases, the signals are returned
#' in the specified order. If \code{length(signals) == 0} or \code{is.na(signals)}, then no records are read
#' and only the file header is returned.
#' @param multifreq Character string. Specifies what to do when signals with different sampling frequencies are
#' present in the file. Options are:
#' \describe{
#' \item{\code{"group"}}{Signals with identical sampling frequencies are grouped together in a single \code{ctd} object
#' (Default)}
#' \item{\code{"separate"}}{Each signal will have its own \code{ctd} object}
#' \item{\code{"downsample"}}{Signals with higher sampling frequencies are downsampled to the lowest frequency in the file}
#' \item{\code{"upsample"}}{Signals with lower sampling frequencies are upsampled to the highest frequency in the file}
#' }
#' Up- and downsampling are performed by the \code{resample} function in the \code{\link{gsignal}} package. Note that this can
#' take a lot of time, especially in case of upsampling to the frequency of the annotation signal (usually 1-2 kHz).
#' @param physical Logical. If true, convert the digital data values to physical values. \eqn{physval = gain * digival + offset},
#' where \eqn{gain = (physmax - physmin) / (digimax - digimin)}, and \eqn{offset = physmax - gain * digimax}. See the description
#' of the header below for the definition of the variables used in these equations.
#'
#' @return \code{list (head, ctd, ..., anno)}. The function returns a \code{list} containing three possible elements: the
#' file header, the data, and annotations, if present. The header is always returned; the other two elements depend on the
#' input arguments.
#' \describe{
#' \item{\strong{header}}{A list named \code{'header'} consisting of the following elements:
#' \describe{
#' \item{General header:\cr}{
#' \describe{
#' \item{\code{id0}}{First identification byte. For BDF: -1 (0xFF), for EDF: ASCII "0" (0x30, old) or " " (0x32, new)}
#' \item{\code{id1 }}{Identification string, "" for EDF, "BIOSEMI" for BDF (7 bytes ASCII)}
#' \item{\code{lsi}}{Local subject identification (80 bytes ASCII)}
#' \item{\code{lri}}{Local recording identification (80 bytes ASCII)}
#' \item{\code{date}}{Starting date of recording (8 bytes ASCII - "dd.mm.yy")}
#' \item{\code{time}}{Starting time of recording (8 bytes ASCII - "hh.mm.ss")}
#' \item{\code{nbytes}}{Number of bytes in header record (converted from 8 bytes ASCII to numeric)}
#' \item{\code{version}}{Reserved. Used for version of format. Empty for EDF, "EDF+C" for EDF continuous data,
#' "EDF+D" for EDF discontinuous data, "BIOSEMI" or "24BIT" for BDF (44 bytes ASCII)}
#' \item{\code{nrec}}{Number of data records (converted from 8 bytes ASCII to numeric), -1 if unknown}
#' \item{\code{duration}}{Duration of one data record in seconds (converted from 8 bytes ASCII to numeric)}
#' \item{\code{ns}}{Number of signals in each data record (converted from 4 bytes ASCII to numeric)}
#' }
#' }
#' \item{Signal-specific header (one for each signal in the file):\cr}{
#' \describe{
#' \item{\code{label}}{Signal labels, e.g. "Fpz", "Fz" (ns * 16 bytes ASCII)}
#' \item{\code{transducer}}{Transducer type, e.g. "AgAgCl electrode" (ns * 80 bytes ASCII)}
#' \item{\code{dimension}}{Physical dimension of signal, e.g., "uV", "Ohm" (ns * 8 bytes ASCII)}
#' \item{\code{physmin}}{Physical minimum in units of dimension (ns * 8 bytes ASCII converted to numeric)}
#' \item{\code{physmax}}{Physical maximum in units of dimension (ns * 8 bytes ASCII converted to numeric)}
#' \item{\code{digimin}}{Digital minimum ((ns * 8 bytes ASCII converted to numeric)}
#' \item{\code{digimax}}{Digital maximum (ns * 8 bytes ASCII converted to numeric)}
#' \item{\code{gain}}{Gain factor used for mapping digital to physical values (ns * numeric)}
#' \item{\code{offset}}{Offset value used for mapping digital to physical values (ns * numeric)}
#' \item{\code{prefilt}}{Prefiltering, e.g., "TC:3s;LP:70Hz" (ns * 80 bytes ASCII)}
#' \item{\code{nsamp}}{Number of samples in one data record (ns * 8 bytes ASCII converted to numeric)}
#' \item{\code{resvd}}{Reserved (ns * 32 bytes ASCII)}
#' }
#' }
#' }
#' }
#' \item{\strong{data}}{One or more Continuous Time Domain (\code{\link{ctd}}) objects. If the input file contains signals
#' which all have the same sampling frequency, or if either the \code{"upsample"} or \code{"downsample"} options were
#' specified in the \code{multifreq} parameter, then there is a single \code{ctd} object, and it is named '\code{ctd}'.\cr\cr
#' In case of multiple sampling frequencies, then if the \code{multifreq = "group"}, as many \code{ctd} objects are
#' returned as there are sampling frequencies, and they are named \code{ctdF}, where F equals the sampling frequency
#' of the data in that object. For instance, for a data file with sampling frequencies of 256 and 1024 Hz, two \code{ctd}
#' objects will be returned, named \code{ctd256} and \code{ctd1024}.\cr\cr
#' In case \code{multifreq = "separate"} is specified, there will be as many \code{ctd} objects as there are signals,
#' and they will be named \code{'ctd1...ctdN'}, where \code{N} equals the number of signals in the file.
#' }
#' \item{\strong{annotations}}{A (list of) data frame(s) named \code{'anno'} containing the annotations parsed from
#' the EDF+/BDF+ annotation channels, or from the Biosemi 'status' channel. Each annotation consists of the following fields:
#' \describe{
#' \item{\code{sample}}{sample number of annotation}
#' \item{\code{time}}{corresponding time, relative to beginning of file}
#' \item{\code{value}}{value of annotation, as per EDF+ specification, see \url{https://www.edfplus.info/specs/edfplus.html}.
#' For Biosemi this ia a 16 bit number, as only bits 0-15 of 24 are used as triggers, see
#' \url{http://biosemi.com/faq/trigger_signals.htm}}.
#' }
#' For EDF+ and BDF+ files, the Time-stamped Annotation Lists (TALs) will also be returned. In this case, sample numbers and
#' times are relative to the beginning of the file, irrespective of the records read.
#' }
#' }
#'
#' @examples
#' \dontrun{
#' edf <- readEDF('test.edf', signals = 1:10, records = "all")
#' summary(edf$ctd)
#' header <- readEDF(file.bdf, records = NA)
#' header
#' }
#' @importFrom mmap mmap munmap int16 int24 char
#' @importFrom gsignal resample
#' @export
# file <- "/media/geert/DATA/Data/Stefan/Robin/edf/p05.edf"
# signals=28:30
# signals<-c("F3","Fz","F4")
# #signals=c(1:21, 23,24)
# #signals=22
# records=NULL
readEDF <- function (file, records = "all", signals = "all",
multifreq = c("group", "separate", "downsample", "upsample"),
physical = TRUE) {
# helper function to trim leading and trailing spaces from header strings
trim.spaces <- function(a) gsub("^ *","",gsub(" *$","",a))
# check parameters
if (missing(file) || is.null(file) || !is.character(file)) stop(paste("Invalid 'file' parameter:", file))
if (!file.exists(file)) stop(paste(file, "does not exist"))
# Open input file and read header variables
fp <- file(file, "rb")
if (!isOpen(fp)) stop(paste("Unable to open input file:", file))
result <- tryCatch({
id0 <- readBin(fp, integer(), n = 1L, size = 1L, endian = "little")
id1 <- readChar(fp, 7)
lsi <- trim.spaces(readChar(fp, 80))
lri <- trim.spaces(readChar(fp, 80))
date <- readChar(fp, 8)
time <- readChar(fp, 8)
nbytes <- as.numeric(readChar(fp, 8))
version <- trim.spaces(readChar(fp, 44))
nrec <- as.numeric(readChar(fp, 8))
if (nrec <= 0) stop(paste("\nInvalid number of records, nrec =", nrec))
duration <- as.numeric(readChar(fp, 8))
ns <- as.numeric(readChar(fp, 4))
if (ns <= 0) stop(paste("\nInvalid number of signals, ns =", ns))
label <- array('', ns)
for (is in seq_len(ns)) {
label[is] <- trim.spaces(readChar(fp, 16))
}
transducer <- array('', ns)
for (is in seq_len(ns)) {
transducer[is] <- trim.spaces(readChar(fp, 80))
}
dimension <- array('', ns)
for (is in seq_len(ns)) {
dimension[is] <- trim.spaces(readChar(fp, 8))
}
physmin <- array(0, ns)
for (is in seq_len(ns)) {
physmin[is] <- as.numeric(readChar(fp, 8))
}
physmax <- array(0, ns)
for (is in seq_len(ns)) {
physmax[is] <- as.numeric(readChar(fp, 8))
}
digimin <- array(0, ns)
for (is in seq_len(ns)) {
digimin[is] <- as.numeric(readChar(fp, 8))
}
digimax <- array(0, ns)
for (is in seq_len(ns)) {
digimax[is] <- as.numeric(readChar(fp, 8))
}
gain <- (physmax - physmin) / (digimax - digimin)
offset <- physmax - gain * digimax
prefilter <- array('', ns)
for (is in seq_len(ns)) {
prefilter[is] <- trim.spaces(readChar(fp, 80))
}
nsamp <- array(0, ns)
for (is in seq_len(ns)) {
nsamp[is] <- as.numeric(readChar(fp, 8))
}
resvd <- array(0, ns)
for (is in seq_len(ns)) {
resvd[is] <- trim.spaces(readChar(fp, 32))
}
}, #end of tryCatch expression
warning = function(w) warning(paste("Warning reading header:", w)),
error = function(e) stop(paste("Error reading header:", e))
) #end of tryCatch block
# check header size
if (nbytes != ns*256 + 256) warning("inconsistent header size, attempting to continue...")
head <- list(id0 = id0, id1 = id1, lsi = lsi, lri = lri, date = date, time = time,
nbytes = nbytes, version = version, nrec = nrec, duration = duration,
ns = ns, label = label, transducer = transducer, dimension = dimension,
physmin = physmin, physmax = physmax, digimin = digimin,
digimax = digimax, gain = gain, offset = offset, prefilter = prefilter,
nsamp = nsamp, resvd = resvd)
# Determine how many records and signals to read. Return header if either
# no records or no signals are requested
if (missing(records) || is.null(records) || trim.spaces(records) == 'all') records <- seq_len(nrec)
if (missing(signals) || is.null(signals) || trim.spaces(signals) == 'all') signals <- seq_len(ns)
if (length(records) == 0 || (length(records) == 1 && (is.na(records) || records == 0)) ||
length(signals) == 0 || (length(signals) == 1 && (is.na(signals) || signals == 0))) {
close(fp)
return(head)
}
# Translate named signals as numbers
if (typeof(signals) == "character") signals <- match(signals, label)
# Avoid array out of bounds conditions
records <- unique(abs(records[!is.na(records) & records > 0 & records <= nrec]))
signals <- unique(abs(signals[!is.na(signals) & signals > 0 & signals <= ns]))
# Read all samples as raw() and write them back to a temporary file which will
# then be mapped into memory by the mmap function. This is a bit clumsy,
# but using mmap is the fastest way to read 24-bit integers. Unfortunately,
# I have found no way to map the entire original file (including the header) into memory.
# The mmap function does have an offset variable, but it can only be equal to the system pagesize.
# So, the solution is to write back the samples to disk withoug the header, and to map this temporary
# file to memory starting with the first byte. Tests have shown that this is faster
# than fooling around with bytes in a big loop.
# Determine whether to read 16-bit (EDF) or 24-bit (BDF) integers
size <- ifelse(id0==-1 && substr(id1,1,7)=='BIOSEMI' &&
(substr(version,1,5)=='24BIT' || substr(version,1,3) == 'BDF'), 3, 2)
# open temporary file for writing single bytes
tmpf <- tempfile()
tmp <- file(tmpf, "wb")
if (!isOpen(tmp)) {
message (paste("Unable to open temporary file for writing:", tmpf))
close (fp)
return(head)
}
# read everything at once
bytes2read <- size * sum(nsamp) * nrec
sample <- readBin (fp, raw(), n=bytes2read, size=1L, endian="little")
writeBin (sample, tmp, raw())
close (fp)
close (tmp)
# Map the temporary file into memory, either as 24-bit int or as 16-bit int
if (size == 3) {
smp <- mmap::mmap(tmpf, mode=mmap::int24())
} else {
smp <- mmap::mmap(tmpf, mode=mmap::int16())
}
recsize <- sum(nsamp) # total record size including the signals we do not want to read
if (length(smp) != recsize * nrec) { # check against length of smp
message (paste("Header/data mismatch",
"\nHeader indicates to read", recsize * nrec, "bytes,",
"\nbut", length(smp), "bytes were read"))
return(head)
}
# Get the signals we want from the records we want and move them into temporay arrays
# these arrays can be of different length
dat <- ptr <- list()
for (s in signals) {
dat[[s]] <- array(0, length(records) * nsamp[s])
ptr[[s]] <- 1
}
for (r in records) {
for (s in signals) {
smp.f <- 1 + (r - 1) * recsize + ifelse(s == 1, 0, sum(nsamp[1:(s-1)]))
smp.t <- smp.f + nsamp[s] - 1
dat.f <- ptr[[s]]
dat.t <- ptr[[s]] + nsamp[s] - 1
dat[[s]][dat.f:dat.t] <- smp[smp.f:smp.t]
ptr[[s]] <- ptr[[s]] + nsamp[s]
}
}
# convert digital to physical values, if requested
if (physical) {
for (s in signals) {
if (label[s] != 'EDF Annotations' && label[s] != 'BDF Annotations' && label[s] != 'Status') {
dat[[s]] <- dat[[s]] * gain[s] + offset[s]
}
}
}
# Determine how the data part of the output should look like.
multifreq <- match.arg(multifreq)
# 1. If there is a single sampling frequency, then make a single ctd object
freq <- nsamp / duration
frex <- unique(freq[signals])
nfrex <- length(frex)
if (nfrex == 1) {
ctd <- ctd(matrix(unlist(dat[signals]), nrow = frex * duration * length(records), ncol = length(signals)), frex)
colnames(ctd)[1:ns(ctd)] <- label[signals]
}
# 2. If there are multiple frequencies and multifreq == 'group', then group all signals with
# the same sampling frequencies into one ctd object. Name the ctd objects ctdF where F is the
# sampling frequency.
if (nfrex > 1 && multifreq == 'group') {
ctd <- mtx <- list()
m <- match(freq, frex)[signals] #match the selected signals to the frex
tab <- table(m) #number of signal per freq
ptr <- array(1, nfrex) # keep track of column
for (f in seq_along(frex)) {
mtx[[f]] <- matrix(0, nrow = frex[f] * duration * length(records), ncol = tab[f])
cn <- NULL
for (s in signals) {
if (freq[s] == frex[f]) {
mtx[[f]][, ptr[f]] <- dat[[s]]
ptr[f] <- ptr[f] + 1
cn <- c(cn, label[s])
}
}
colnames(mtx[[f]]) <- cn
ctd[[f]] <- ctd(mtx[[f]], frex[f])
}
names(ctd) <- paste0('ctd', frex)
}
# 3. If there are multiple frequencies and multifreq == 'separate', then there will be a ctd object for
# each signal, and they will be named ctd1...ctdN, where N is the number of signals
if (nfrex > 1 && multifreq == 'separate') {
ctd <- list()
ptr <- 1
for (s in signals) {
ctd[[ptr]] <- ctd(dat[[s]], freq[s])
colnames(ctd[[ptr]])[1] <- label[s]
ptr <- ptr + 1
}
names(ctd) <- paste0('ctd', 1:length(ctd))
}
# 4. If there are multiple frequencies and "upsample" is specified, then determine the highest
# sampling rate and upsample every signal to that frequency
# 5. If there are multiple frequencies and "downsample" is specified, then determine the lowest
# sampling rate and downsample every signal to that frequency
if (nfrex > 1 && (multifreq == "upsample" || multifreq == "downsample")) {
newf <- ifelse(multifreq == "upsample", max(frex), min(frex))
mtx <- matrix(0, nrow = newf * duration * length(records), ncol = length(signals))
col <- 1
for (s in signals) {
if (freq[s] != newf) mtx[, col] <- gsignal::resample(as.vector(dat[[s]]), newf, freq[s])
else mtx[, col] <- as.vector(dat[[s]])
col <- col + 1
}
ctd <- ctd(mtx, newf)
colnames(ctd)[1:ns(ctd)] <- label[signals]
}
# Annotations are ALWAYS parsed if present in the file. For EDF files, there is no annotation channel.
# For Biosemi generated BDF files, there is always a 'Status' channel which is the last signal in the file.
#Try to find this channel in the selected signals.
anno <- NULL
if (id0 == -1 && substr(id1, 1, 7) == 'BIOSEMI' && substr(version,1,5)=='24BIT') {
annosig <- which(label[signals] == 'Status')
nanno <- length(annosig)
if (nanno > 1) annosig <- annosig[1]
if (nanno != 0) {
df <- diff(dat[[annosig]])
trig <- which(df>0 & df <= 2^16) # Biosemi uses 16-bit triggers
if (length(trig) > 0) {
as <- trig + 1 # high level is moment of trigger
at <- (trig / (nsamp[annosig] / duration))
av <- df[trig]
anno <- data.frame(sample=as, time=at, value=av);
}
}
} else {
# For B/EDF+C and B/EDF+D files, try to find the annotation channels among the selected signals.
# These signals will contain Time-stamped Annotation Lists (TALs), containing single-byte raw characters.
# Re-open the temporary file as char().
getTALs <- function (chan, sf) {
TALs <- NULL
start0 <- which(chan == 0x2b | chan == 0x2d) #TAL starts with '+' or '-'
# bugfix 20191210
if (length(start0) == 0) return(NULL)
start1 <- start0 - 1 #preceded by 0x00
if (start0[1] == 1) start <- c(1, start0[2:length(start0)][which(chan[start1] == 0)])
else start <- start0[which(chan[start1] == 0)]
end20 <- which(chan == 20) #TAL ends with 20
end0 <- end20 + 1 #followed by 0
end <- end20[which(chan[end0] == 0)] #last sample handled implicitly
n <- length(start)
if (n == length(end)) {
for (tal in seq_along(start)) TALs <- rbind(TALs, rawToChar(chan[start[tal]:end[tal]]))
}
parseTALs(TALs, sf)
}
parseTALs <- function(tals, sf) {
if(nrow(tals) <= 0) return(NULL)
s <- t <- d <- v <- array(0, nrow(tals))
for (r in seq_len(nrow(tals))) {
split1 <- unlist(strsplit(tals[r], split = "\024")) #octal 24 = hex 14 = dec 20
split2 <- unlist(strsplit(split1, split = "\025")) #octal 25 = hex 15 = dec 21
suppressWarnings(t[r] <- as.numeric(split2[1])) #allow NAs by coercion
s[r] <- round(t[r] * sf) + 1
suppressWarnings(d[r] <- as.numeric(split2[2])) #allow NAs by coercion
v[r] <- trim.spaces(split1[2])
}
data.frame(sample = s, time = t, value = v, duration = d, tal = tals,
stringsAsFactors = FALSE)
}
m <- match(label, c('EDF Annotations', 'BDF Annotations'))
annosig <- which(m > 0)
sig <- anno <- list()
if (!is.null(annosig)) {
smp <- mmap::mmap(tmpf, mode=mmap::char())
for (s in seq_along(annosig)) {
sig[[s]] <- raw(length = size * length(records) * nsamp[annosig[s]])
ptr <- 1
for (r in seq_len(nrec)) {
smp.f <- 1 + (r - 1) * recsize * size + sum(nsamp[1:(annosig[s]-1)] * size)
smp.t <- smp.f + nsamp[annosig[s]] * size - 1
dat.f <- ptr
dat.t <- ptr + nsamp[annosig[s]] * size - 1
sig[[s]][dat.f:dat.t] <- smp[smp.f:smp.t]
ptr <- ptr + nsamp[annosig[s]] * size
}
anno[[s]] <- getTALs(sig[[s]], nsamp[annosig[s]] / duration)
}
}
if(length(anno) == 1) anno <- anno[[1]] # no point in returning a list if there is only one anno channel
}
# Clean up
mmap::munmap(smp)
unlink(tmpf)
# Ready
return(list(header = head, ctd = ctd, anno = anno))
}
# s <- sprintf("%04X", ctd$data[1,1])
# h <- sapply(seq(1, nchar(s), by=2), function(x) substr(s, x, x+1))
# rawToChar(as.raw(strtoi(h[1])))
#
# t <- tempfile()
# fp <- file(t, "wb")
# writeBin(ctd$data[1:5,1], fp, numeric())
# close(fp)
# fp <- file(t, "rb")
# asc <- readBin(fp, character(), n=10)
# close(fp)
gjmvanboxtel/eegr documentation built on May 20, 2023, 4:26 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# readEDF.R - Function to read an EDF/BDF file into an R data frame
# Copyright (C) 2012, 2019 Geert van Boxtel,
# Tilburg University, G.J.M.vBoxtel@tilburguniversity.edu
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
# See also: http://www.gnu.org/licenses/gpl-2.0.txt
#
# Version history:
# 20120207 GvB Original version
# 20190701 GvB complete overhaul for eegr v0.1.0
# 20201104 GvB link to gsignal instead of signal package
# 20230516 GvB bugfixes in calculating record size and up/downsampling
#
#-----------------------------------------------------------------------------------------------
#' readEDF
#'
#' Read data from a European Data Format (EDF, EDF+, BDF, BDF+) file.
#'
#' The function readEDF reads files containing electrophysiological data stored in the European Data
#' Format (EDF) or its variants.
#'
#' EDF was published in 1992 and stores multichannel data, allowing different sample rates for each signal.
#' Internally it includes a header and one or more data records. The header contains some general information
#' (patient identification, start time...) and technical specs of each signal (calibration, sampling rate, filtering, ...),
#' coded as ASCII characters. The data records contain samples as little-endian 16-bit integers.
#'
#' EDF+ was published in 2003 and is largely compatible to EDF, but EDF+ files also allow coding discontinuous recordings
#' as well as annotations, stimuli and events in UTF-8 format.
#'
#' Biosemi developed an EEG acquisition system with 24-bit A/D converters, which could not be stored in the 16-bit
#' EDF format. They introduced a 24-bit variant of EDF called BDF (Biosemi Data Format), and a BDF+ variant is sometimes also
#' found.
#'
#' The function \code{readEDF} can read all of these variants and produce different output that can be used for further processing
#' by function in the \code{eegr} package.
#'
#' @references Kemp B., Värri, A., Rosa, A.C., Nielsen, K.D. and Gade, J. (1992). A simple format for exchange of
#' digitized polygraphic recordings. \emph{Electroencephalogr Clin Neurophysiol}, \emph{82(5)}, 391-393.
#' @references Kemp, B. and Olivan, J. (2003). European data format ’plus’ (EDF+), an EDF alike standard format
#' for the exchange of physiological data. \emph{Clin Neurophysiol}, \emph{114(9)}, 1755-1761.
#' @references European Data Format. \url{https://www.edfplus.info/}
#' @references Which file format does Biosemi use? \url{https://www.biosemi.com/faq/file_format.htm}
#'
#' @author Geert van Boxtel, \email{G.J.M.vanBoxtel@@gmail.com}
#'
#' @param file Character string. The name of the file which the data are to be read from. \code{file}
#' can also be a complete URL. (For the supported URL schemes,see the ‘URLs’ section of the help for url.)
#' @param records Numeric or \code{'all'}. In EDF, signals are stored in 'records', usually containing 0.5 or 1
#' second of data. This function either reads all records (default: \code{all}), or a subset of the records, thus
#' limiting the time interval of the signals. \code{records} can be a range (e.g., \code{1:50}) or a random sequence
#' of records (e.g., \code{c(1,5,12,3)}). In the latter case, the records are returned in the specified order (only
#' useful for discontinuous data). If \code{length(records) == 0} or \code{is.na(records)}, then no records are read
#' and only the file header is returned.
#' @param signals Numeric, character string, or \code{'all'}. Specifies the signals to read from the input file.
#' \code{signals} can be a range (e.g., \code{1:21}) a random sequence of signals (e.g., \code{c(1,5,12,3)}),
#' or a sequence of signal labels (e.g., \code{c("F3", "Fz", "F4")}. In the latter two cases, the signals are returned
#' in the specified order. If \code{length(signals) == 0} or \code{is.na(signals)}, then no records are read
#' and only the file header is returned.
#' @param multifreq Character string. Specifies what to do when signals with different sampling frequencies are
#' present in the file. Options are:
#' \describe{
#' \item{\code{"group"}}{Signals with identical sampling frequencies are grouped together in a single \code{ctd} object
#' (Default)}
#' \item{\code{"separate"}}{Each signal will have its own \code{ctd} object}
#' \item{\code{"downsample"}}{Signals with higher sampling frequencies are downsampled to the lowest frequency in the file}
#' \item{\code{"upsample"}}{Signals with lower sampling frequencies are upsampled to the highest frequency in the file}
#' }
#' Up- and downsampling are performed by the \code{resample} function in the \code{\link{gsignal}} package. Note that this can
#' take a lot of time, especially in case of upsampling to the frequency of the annotation signal (usually 1-2 kHz).
#' @param physical Logical. If true, convert the digital data values to physical values. \eqn{physval = gain * digival + offset},
#' where \eqn{gain = (physmax - physmin) / (digimax - digimin)}, and \eqn{offset = physmax - gain * digimax}. See the description
#' of the header below for the definition of the variables used in these equations.
#'
#' @return \code{list (head, ctd, ..., anno)}. The function returns a \code{list} containing three possible elements: the
#' file header, the data, and annotations, if present. The header is always returned; the other two elements depend on the
#' input arguments.
#' \describe{
#' \item{\strong{header}}{A list named \code{'header'} consisting of the following elements:
#' \describe{
#' \item{General header:\cr}{
#' \describe{
#' \item{\code{id0}}{First identification byte. For BDF: -1 (0xFF), for EDF: ASCII "0" (0x30, old) or " " (0x32, new)}
#' \item{\code{id1 }}{Identification string, "" for EDF, "BIOSEMI" for BDF (7 bytes ASCII)}
#' \item{\code{lsi}}{Local subject identification (80 bytes ASCII)}
#' \item{\code{lri}}{Local recording identification (80 bytes ASCII)}
#' \item{\code{date}}{Starting date of recording (8 bytes ASCII - "dd.mm.yy")}
#' \item{\code{time}}{Starting time of recording (8 bytes ASCII - "hh.mm.ss")}
#' \item{\code{nbytes}}{Number of bytes in header record (converted from 8 bytes ASCII to numeric)}
#' \item{\code{version}}{Reserved. Used for version of format. Empty for EDF, "EDF+C" for EDF continuous data,
#' "EDF+D" for EDF discontinuous data, "BIOSEMI" or "24BIT" for BDF (44 bytes ASCII)}
#' \item{\code{nrec}}{Number of data records (converted from 8 bytes ASCII to numeric), -1 if unknown}
#' \item{\code{duration}}{Duration of one data record in seconds (converted from 8 bytes ASCII to numeric)}
#' \item{\code{ns}}{Number of signals in each data record (converted from 4 bytes ASCII to numeric)}
#' }
#' }
#' \item{Signal-specific header (one for each signal in the file):\cr}{
#' \describe{
#' \item{\code{label}}{Signal labels, e.g. "Fpz", "Fz" (ns * 16 bytes ASCII)}
#' \item{\code{transducer}}{Transducer type, e.g. "AgAgCl electrode" (ns * 80 bytes ASCII)}
#' \item{\code{dimension}}{Physical dimension of signal, e.g., "uV", "Ohm" (ns * 8 bytes ASCII)}
#' \item{\code{physmin}}{Physical minimum in units of dimension (ns * 8 bytes ASCII converted to numeric)}
#' \item{\code{physmax}}{Physical maximum in units of dimension (ns * 8 bytes ASCII converted to numeric)}
#' \item{\code{digimin}}{Digital minimum ((ns * 8 bytes ASCII converted to numeric)}
#' \item{\code{digimax}}{Digital maximum (ns * 8 bytes ASCII converted to numeric)}
#' \item{\code{gain}}{Gain factor used for mapping digital to physical values (ns * numeric)}
#' \item{\code{offset}}{Offset value used for mapping digital to physical values (ns * numeric)}
#' \item{\code{prefilt}}{Prefiltering, e.g., "TC:3s;LP:70Hz" (ns * 80 bytes ASCII)}
#' \item{\code{nsamp}}{Number of samples in one data record (ns * 8 bytes ASCII converted to numeric)}
#' \item{\code{resvd}}{Reserved (ns * 32 bytes ASCII)}
#' }
#' }
#' }
#' }
#' \item{\strong{data}}{One or more Continuous Time Domain (\code{\link{ctd}}) objects. If the input file contains signals
#' which all have the same sampling frequency, or if either the \code{"upsample"} or \code{"downsample"} options were
#' specified in the \code{multifreq} parameter, then there is a single \code{ctd} object, and it is named '\code{ctd}'.\cr\cr
#' In case of multiple sampling frequencies, then if the \code{multifreq = "group"}, as many \code{ctd} objects are
#' returned as there are sampling frequencies, and they are named \code{ctdF}, where F equals the sampling frequency
#' of the data in that object. For instance, for a data file with sampling frequencies of 256 and 1024 Hz, two \code{ctd}
#' objects will be returned, named \code{ctd256} and \code{ctd1024}.\cr\cr
#' In case \code{multifreq = "separate"} is specified, there will be as many \code{ctd} objects as there are signals,
#' and they will be named \code{'ctd1...ctdN'}, where \code{N} equals the number of signals in the file.
#' }
#' \item{\strong{annotations}}{A (list of) data frame(s) named \code{'anno'} containing the annotations parsed from
#' the EDF+/BDF+ annotation channels, or from the Biosemi 'status' channel. Each annotation consists of the following fields:
#' \describe{
#' \item{\code{sample}}{sample number of annotation}
#' \item{\code{time}}{corresponding time, relative to beginning of file}
#' \item{\code{value}}{value of annotation, as per EDF+ specification, see \url{https://www.edfplus.info/specs/edfplus.html}.
#' For Biosemi this ia a 16 bit number, as only bits 0-15 of 24 are used as triggers, see
#' \url{http://biosemi.com/faq/trigger_signals.htm}}.
#' }
#' For EDF+ and BDF+ files, the Time-stamped Annotation Lists (TALs) will also be returned. In this case, sample numbers and
#' times are relative to the beginning of the file, irrespective of the records read.
#' }
#' }
#'
#' @examples
#' \dontrun{
#' edf <- readEDF('test.edf', signals = 1:10, records = "all")
#' summary(edf$ctd)
#' header <- readEDF(file.bdf, records = NA)
#' header
#' }
#' @importFrom mmap mmap munmap int16 int24 char
#' @importFrom gsignal resample
#' @export
# file <- "/media/geert/DATA/Data/Stefan/Robin/edf/p05.edf"
# signals=28:30
# signals<-c("F3","Fz","F4")
# #signals=c(1:21, 23,24)
# #signals=22
# records=NULL
readEDF <- function (file, records = "all", signals = "all",
multifreq = c("group", "separate", "downsample", "upsample"),
physical = TRUE) {
# helper function to trim leading and trailing spaces from header strings
trim.spaces <- function(a) gsub("^ *","",gsub(" *$","",a))
# check parameters
if (missing(file) || is.null(file) || !is.character(file)) stop(paste("Invalid 'file' parameter:", file))
if (!file.exists(file)) stop(paste(file, "does not exist"))
# Open input file and read header variables
fp <- file(file, "rb")
if (!isOpen(fp)) stop(paste("Unable to open input file:", file))
result <- tryCatch({
id0 <- readBin(fp, integer(), n = 1L, size = 1L, endian = "little")
id1 <- readChar(fp, 7)
lsi <- trim.spaces(readChar(fp, 80))
lri <- trim.spaces(readChar(fp, 80))
date <- readChar(fp, 8)
time <- readChar(fp, 8)
nbytes <- as.numeric(readChar(fp, 8))
version <- trim.spaces(readChar(fp, 44))
nrec <- as.numeric(readChar(fp, 8))
if (nrec <= 0) stop(paste("\nInvalid number of records, nrec =", nrec))
duration <- as.numeric(readChar(fp, 8))
ns <- as.numeric(readChar(fp, 4))
if (ns <= 0) stop(paste("\nInvalid number of signals, ns =", ns))
label <- array('', ns)
for (is in seq_len(ns)) {
label[is] <- trim.spaces(readChar(fp, 16))
}
transducer <- array('', ns)
for (is in seq_len(ns)) {
transducer[is] <- trim.spaces(readChar(fp, 80))
}
dimension <- array('', ns)
for (is in seq_len(ns)) {
dimension[is] <- trim.spaces(readChar(fp, 8))
}
physmin <- array(0, ns)
for (is in seq_len(ns)) {
physmin[is] <- as.numeric(readChar(fp, 8))
}
physmax <- array(0, ns)
for (is in seq_len(ns)) {
physmax[is] <- as.numeric(readChar(fp, 8))
}
digimin <- array(0, ns)
for (is in seq_len(ns)) {
digimin[is] <- as.numeric(readChar(fp, 8))
}
digimax <- array(0, ns)
for (is in seq_len(ns)) {
digimax[is] <- as.numeric(readChar(fp, 8))
}
gain <- (physmax - physmin) / (digimax - digimin)
offset <- physmax - gain * digimax
prefilter <- array('', ns)
for (is in seq_len(ns)) {
prefilter[is] <- trim.spaces(readChar(fp, 80))
}
nsamp <- array(0, ns)
for (is in seq_len(ns)) {
nsamp[is] <- as.numeric(readChar(fp, 8))
}
resvd <- array(0, ns)
for (is in seq_len(ns)) {
resvd[is] <- trim.spaces(readChar(fp, 32))
}
}, #end of tryCatch expression
warning = function(w) warning(paste("Warning reading header:", w)),
error = function(e) stop(paste("Error reading header:", e))
) #end of tryCatch block
# check header size
if (nbytes != ns*256 + 256) warning("inconsistent header size, attempting to continue...")
head <- list(id0 = id0, id1 = id1, lsi = lsi, lri = lri, date = date, time = time,
nbytes = nbytes, version = version, nrec = nrec, duration = duration,
ns = ns, label = label, transducer = transducer, dimension = dimension,
physmin = physmin, physmax = physmax, digimin = digimin,
digimax = digimax, gain = gain, offset = offset, prefilter = prefilter,
nsamp = nsamp, resvd = resvd)
# Determine how many records and signals to read. Return header if either
# no records or no signals are requested
if (missing(records) || is.null(records) || trim.spaces(records) == 'all') records <- seq_len(nrec)
if (missing(signals) || is.null(signals) || trim.spaces(signals) == 'all') signals <- seq_len(ns)
if (length(records) == 0 || (length(records) == 1 && (is.na(records) || records == 0)) ||
length(signals) == 0 || (length(signals) == 1 && (is.na(signals) || signals == 0))) {
close(fp)
return(head)
}
# Translate named signals as numbers
if (typeof(signals) == "character") signals <- match(signals, label)
# Avoid array out of bounds conditions
records <- unique(abs(records[!is.na(records) & records > 0 & records <= nrec]))
signals <- unique(abs(signals[!is.na(signals) & signals > 0 & signals <= ns]))
# Read all samples as raw() and write them back to a temporary file which will
# then be mapped into memory by the mmap function. This is a bit clumsy,
# but using mmap is the fastest way to read 24-bit integers. Unfortunately,
# I have found no way to map the entire original file (including the header) into memory.
# The mmap function does have an offset variable, but it can only be equal to the system pagesize.
# So, the solution is to write back the samples to disk withoug the header, and to map this temporary
# file to memory starting with the first byte. Tests have shown that this is faster
# than fooling around with bytes in a big loop.
# Determine whether to read 16-bit (EDF) or 24-bit (BDF) integers
size <- ifelse(id0==-1 && substr(id1,1,7)=='BIOSEMI' &&
(substr(version,1,5)=='24BIT' || substr(version,1,3) == 'BDF'), 3, 2)
# open temporary file for writing single bytes
tmpf <- tempfile()
tmp <- file(tmpf, "wb")
if (!isOpen(tmp)) {
message (paste("Unable to open temporary file for writing:", tmpf))
close (fp)
return(head)
}
# read everything at once
bytes2read <- size * sum(nsamp) * nrec
sample <- readBin (fp, raw(), n=bytes2read, size=1L, endian="little")
writeBin (sample, tmp, raw())
close (fp)
close (tmp)
# Map the temporary file into memory, either as 24-bit int or as 16-bit int
if (size == 3) {
smp <- mmap::mmap(tmpf, mode=mmap::int24())
} else {
smp <- mmap::mmap(tmpf, mode=mmap::int16())
}
recsize <- sum(nsamp) # total record size including the signals we do not want to read
if (length(smp) != recsize * nrec) { # check against length of smp
message (paste("Header/data mismatch",
"\nHeader indicates to read", recsize * nrec, "bytes,",
"\nbut", length(smp), "bytes were read"))
return(head)
}
# Get the signals we want from the records we want and move them into temporay arrays
# these arrays can be of different length
dat <- ptr <- list()
for (s in signals) {
dat[[s]] <- array(0, length(records) * nsamp[s])
ptr[[s]] <- 1
}
for (r in records) {
for (s in signals) {
smp.f <- 1 + (r - 1) * recsize + ifelse(s == 1, 0, sum(nsamp[1:(s-1)]))
smp.t <- smp.f + nsamp[s] - 1
dat.f <- ptr[[s]]
dat.t <- ptr[[s]] + nsamp[s] - 1
dat[[s]][dat.f:dat.t] <- smp[smp.f:smp.t]
ptr[[s]] <- ptr[[s]] + nsamp[s]
}
}
# convert digital to physical values, if requested
if (physical) {
for (s in signals) {
if (label[s] != 'EDF Annotations' && label[s] != 'BDF Annotations' && label[s] != 'Status') {
dat[[s]] <- dat[[s]] * gain[s] + offset[s]
}
}
}
# Determine how the data part of the output should look like.
multifreq <- match.arg(multifreq)
# 1. If there is a single sampling frequency, then make a single ctd object
freq <- nsamp / duration
frex <- unique(freq[signals])
nfrex <- length(frex)
if (nfrex == 1) {
ctd <- ctd(matrix(unlist(dat[signals]), nrow = frex * duration * length(records), ncol = length(signals)), frex)
colnames(ctd)[1:ns(ctd)] <- label[signals]
}
# 2. If there are multiple frequencies and multifreq == 'group', then group all signals with
# the same sampling frequencies into one ctd object. Name the ctd objects ctdF where F is the
# sampling frequency.
if (nfrex > 1 && multifreq == 'group') {
ctd <- mtx <- list()
m <- match(freq, frex)[signals] #match the selected signals to the frex
tab <- table(m) #number of signal per freq
ptr <- array(1, nfrex) # keep track of column
for (f in seq_along(frex)) {
mtx[[f]] <- matrix(0, nrow = frex[f] * duration * length(records), ncol = tab[f])
cn <- NULL
for (s in signals) {
if (freq[s] == frex[f]) {
mtx[[f]][, ptr[f]] <- dat[[s]]
ptr[f] <- ptr[f] + 1
cn <- c(cn, label[s])
}
}
colnames(mtx[[f]]) <- cn
ctd[[f]] <- ctd(mtx[[f]], frex[f])
}
names(ctd) <- paste0('ctd', frex)
}
# 3. If there are multiple frequencies and multifreq == 'separate', then there will be a ctd object for
# each signal, and they will be named ctd1...ctdN, where N is the number of signals
if (nfrex > 1 && multifreq == 'separate') {
ctd <- list()
ptr <- 1
for (s in signals) {
ctd[[ptr]] <- ctd(dat[[s]], freq[s])
colnames(ctd[[ptr]])[1] <- label[s]
ptr <- ptr + 1
}
names(ctd) <- paste0('ctd', 1:length(ctd))
}
# 4. If there are multiple frequencies and "upsample" is specified, then determine the highest
# sampling rate and upsample every signal to that frequency
# 5. If there are multiple frequencies and "downsample" is specified, then determine the lowest
# sampling rate and downsample every signal to that frequency
if (nfrex > 1 && (multifreq == "upsample" || multifreq == "downsample")) {
newf <- ifelse(multifreq == "upsample", max(frex), min(frex))
mtx <- matrix(0, nrow = newf * duration * length(records), ncol = length(signals))
col <- 1
for (s in signals) {
if (freq[s] != newf) mtx[, col] <- gsignal::resample(as.vector(dat[[s]]), newf, freq[s])
else mtx[, col] <- as.vector(dat[[s]])
col <- col + 1
}
ctd <- ctd(mtx, newf)
colnames(ctd)[1:ns(ctd)] <- label[signals]
}
# Annotations are ALWAYS parsed if present in the file. For EDF files, there is no annotation channel.
# For Biosemi generated BDF files, there is always a 'Status' channel which is the last signal in the file.
#Try to find this channel in the selected signals.
anno <- NULL
if (id0 == -1 && substr(id1, 1, 7) == 'BIOSEMI' && substr(version,1,5)=='24BIT') {
annosig <- which(label[signals] == 'Status')
nanno <- length(annosig)
if (nanno > 1) annosig <- annosig[1]
if (nanno != 0) {
df <- diff(dat[[annosig]])
trig <- which(df>0 & df <= 2^16) # Biosemi uses 16-bit triggers
if (length(trig) > 0) {
as <- trig + 1 # high level is moment of trigger
at <- (trig / (nsamp[annosig] / duration))
av <- df[trig]
anno <- data.frame(sample=as, time=at, value=av);
}
}
} else {
# For B/EDF+C and B/EDF+D files, try to find the annotation channels among the selected signals.
# These signals will contain Time-stamped Annotation Lists (TALs), containing single-byte raw characters.
# Re-open the temporary file as char().
getTALs <- function (chan, sf) {
TALs <- NULL
start0 <- which(chan == 0x2b | chan == 0x2d) #TAL starts with '+' or '-'
# bugfix 20191210
if (length(start0) == 0) return(NULL)
start1 <- start0 - 1 #preceded by 0x00
if (start0[1] == 1) start <- c(1, start0[2:length(start0)][which(chan[start1] == 0)])
else start <- start0[which(chan[start1] == 0)]
end20 <- which(chan == 20) #TAL ends with 20
end0 <- end20 + 1 #followed by 0
end <- end20[which(chan[end0] == 0)] #last sample handled implicitly
n <- length(start)
if (n == length(end)) {
for (tal in seq_along(start)) TALs <- rbind(TALs, rawToChar(chan[start[tal]:end[tal]]))
}
parseTALs(TALs, sf)
}
parseTALs <- function(tals, sf) {
if(nrow(tals) <= 0) return(NULL)
s <- t <- d <- v <- array(0, nrow(tals))
for (r in seq_len(nrow(tals))) {
split1 <- unlist(strsplit(tals[r], split = "\024")) #octal 24 = hex 14 = dec 20
split2 <- unlist(strsplit(split1, split = "\025")) #octal 25 = hex 15 = dec 21
suppressWarnings(t[r] <- as.numeric(split2[1])) #allow NAs by coercion
s[r] <- round(t[r] * sf) + 1
suppressWarnings(d[r] <- as.numeric(split2[2])) #allow NAs by coercion
v[r] <- trim.spaces(split1[2])
}
data.frame(sample = s, time = t, value = v, duration = d, tal = tals,
stringsAsFactors = FALSE)
}
m <- match(label, c('EDF Annotations', 'BDF Annotations'))
annosig <- which(m > 0)
sig <- anno <- list()
if (!is.null(annosig)) {
smp <- mmap::mmap(tmpf, mode=mmap::char())
for (s in seq_along(annosig)) {
sig[[s]] <- raw(length = size * length(records) * nsamp[annosig[s]])
ptr <- 1
for (r in seq_len(nrec)) {
smp.f <- 1 + (r - 1) * recsize * size + sum(nsamp[1:(annosig[s]-1)] * size)
smp.t <- smp.f + nsamp[annosig[s]] * size - 1
dat.f <- ptr
dat.t <- ptr + nsamp[annosig[s]] * size - 1
sig[[s]][dat.f:dat.t] <- smp[smp.f:smp.t]
ptr <- ptr + nsamp[annosig[s]] * size
}
anno[[s]] <- getTALs(sig[[s]], nsamp[annosig[s]] / duration)
}
}
if(length(anno) == 1) anno <- anno[[1]] # no point in returning a list if there is only one anno channel
}
# Clean up
mmap::munmap(smp)
unlink(tmpf)
# Ready
return(list(header = head, ctd = ctd, anno = anno))
}
# s <- sprintf("%04X", ctd$data[1,1])
# h <- sapply(seq(1, nchar(s), by=2), function(x) substr(s, x, x+1))
# rawToChar(as.raw(strtoi(h[1])))
#
# t <- tempfile()
# fp <- file(t, "wb")
# writeBin(ctd$data[1:5,1], fp, numeric())
# close(fp)
# fp <- file(t, "rb")
# asc <- readBin(fp, character(), n=10)
# close(fp)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.