Title:  Signal and Image Processing Toolbox for Analyzing Intracranial Electroencephalography Data 

Description:  Implemented fast and memoryefficient Notchfilter, Welchperiodogram, discrete wavelet spectrogram for minutes of highresolution signals, fast 3D convolution, image registration, 3D mesh manipulation; providing fundamental toolbox for intracranial Electroencephalography (iEEG) pipelines. Documentation and examples about 'RAVE' project are provided at <https://openwetware.org/wiki/RAVE>, and the paper by John F. Magnotti, Zhengjia Wang, Michael S. Beauchamp (2020) <doi:10.1016/j.neuroimage.2020.117341>; see 'citation("ravetools")' for details. 
Authors:  Zhengjia Wang [aut, cre], John Magnotti [aut], Michael Beauchamp [aut], Trustees of the University of Pennsylvania [cph] (All files in this package unless explicitly stated in the file or listed in the 'Copyright' section below.), Karim Rahim [cph, ctb] (Contributed to src/ffts.h and stc/ffts.cpp), Thomas Possidente [cph, ctb] (Contributed to R/multitaper.R), Michael Prerau [cph, ctb] (Contributed to R/multitaper.R), Marcus Geelnard [ctb, cph] (TinyThread library, tinythreadpp.bitsnbites.eu, located at inst/include/tthread/), Stefan Schlager [ctb, cph] (Rvcg interface, located at src/vcgCommon.h), Visual Computing Lab, ISTI [ctb, cph] (Copyright holder of vcglib, located at src/vcglib/) 
Maintainer:  Zhengjia Wang <[email protected]> 
License:  GPL (>= 2) 
Version:  0.1.8 
Built:  20240903 20:16:23 UTC 
Source:  https://github.com/dipterix/ravetools 
Bandpass signals
band_pass1(x, sample_rate, lb, ub, domain = 1, ...) band_pass2( x, sample_rate, lb, ub, order, method = c("fir", "butter"), direction = c("both", "forward", "backward"), window = "hamming", ... )
band_pass1(x, sample_rate, lb, ub, domain = 1, ...) band_pass2( x, sample_rate, lb, ub, order, method = c("fir", "butter"), direction = c("both", "forward", "backward"), window = "hamming", ... )
x 
input signals, numeric vector or matrix. 
sample_rate 
sampling frequency 
lb 
lower frequency bound of the bandpassing filter, must be positive 
ub 
upper frequency bound of the bandpassing filter, must be greater than the lower bound and smaller than the half of sampling frequency 
domain 
1 if 
... 
ignored 
order 
the order of the filter, must be positive integer and be less than onethird of the sample rate 
method 
filter type, choices are 
direction 
filter direction, choices are 
window 
window type, can be a character, a function, or a vector.
For character, 
Filtered signals, vector if x
is a vector, or matrix of
the same dimension as x
t < seq(0, 1, by = 0.0005) x < sin(t * 0.4 * pi) + sin(t * 4 * pi) + 2 * sin(t * 120 * pi) oldpar < par(mfrow = c(2, 2), mar = c(3.1, 2.1, 3.1, 0.1)) #  Using band_pass1  y1 < band_pass1(x, 2000, 0.1, 1) y2 < band_pass1(x, 2000, 1, 5) y3 < band_pass1(x, 2000, 10, 80) plot(t, x, type = 'l', xlab = "Time", ylab = "", main = "Mixture of 0.2, 2, and 60Hz") lines(t, y1, col = 'red') lines(t, y2, col = 'blue') lines(t, y3, col = 'green') legend( "topleft", c("Input", "Pass: 0.11Hz", "Pass 15Hz", "Pass 1080Hz"), col = c(par("fg"), "red", "blue", "green"), lty = 1, cex = 0.6 ) # plot pwelch pwelch(x, fs = 2000, window = 4000, noverlap = 2000, plot = 1) pwelch(y1, fs = 2000, window = 4000, noverlap = 2000, plot = 2, col = "red") pwelch(y2, fs = 2000, window = 4000, noverlap = 2000, plot = 2, col = "blue") pwelch(y3, fs = 2000, window = 4000, noverlap = 2000, plot = 2, col = "green") #  Using band_pass2 with FIR filters  order < floor(2000 / 3) z1 < band_pass2(x, 2000, 0.1, 1, method = "fir", order = order) z2 < band_pass2(x, 2000, 1, 5, method = "fir", order = order) z3 < band_pass2(x, 2000, 10, 80, method = "fir", order = order) plot(t, x, type = 'l', xlab = "Time", ylab = "", main = "Mixture of 0.2, 2, and 60Hz") lines(t, z1, col = 'red') lines(t, z2, col = 'blue') lines(t, z3, col = 'green') legend( "topleft", c("Input", "Pass: 0.11Hz", "Pass 15Hz", "Pass 1080Hz"), col = c(par("fg"), "red", "blue", "green"), lty = 1, cex = 0.6 ) # plot pwelch pwelch(x, fs = 2000, window = 4000, noverlap = 2000, plot = 1) pwelch(z1, fs = 2000, window = 4000, noverlap = 2000, plot = 2, col = "red") pwelch(z2, fs = 2000, window = 4000, noverlap = 2000, plot = 2, col = "blue") pwelch(z3, fs = 2000, window = 4000, noverlap = 2000, plot = 2, col = "green") #  Clean this demo  par(oldpar)
t < seq(0, 1, by = 0.0005) x < sin(t * 0.4 * pi) + sin(t * 4 * pi) + 2 * sin(t * 120 * pi) oldpar < par(mfrow = c(2, 2), mar = c(3.1, 2.1, 3.1, 0.1)) #  Using band_pass1  y1 < band_pass1(x, 2000, 0.1, 1) y2 < band_pass1(x, 2000, 1, 5) y3 < band_pass1(x, 2000, 10, 80) plot(t, x, type = 'l', xlab = "Time", ylab = "", main = "Mixture of 0.2, 2, and 60Hz") lines(t, y1, col = 'red') lines(t, y2, col = 'blue') lines(t, y3, col = 'green') legend( "topleft", c("Input", "Pass: 0.11Hz", "Pass 15Hz", "Pass 1080Hz"), col = c(par("fg"), "red", "blue", "green"), lty = 1, cex = 0.6 ) # plot pwelch pwelch(x, fs = 2000, window = 4000, noverlap = 2000, plot = 1) pwelch(y1, fs = 2000, window = 4000, noverlap = 2000, plot = 2, col = "red") pwelch(y2, fs = 2000, window = 4000, noverlap = 2000, plot = 2, col = "blue") pwelch(y3, fs = 2000, window = 4000, noverlap = 2000, plot = 2, col = "green") #  Using band_pass2 with FIR filters  order < floor(2000 / 3) z1 < band_pass2(x, 2000, 0.1, 1, method = "fir", order = order) z2 < band_pass2(x, 2000, 1, 5, method = "fir", order = order) z3 < band_pass2(x, 2000, 10, 80, method = "fir", order = order) plot(t, x, type = 'l', xlab = "Time", ylab = "", main = "Mixture of 0.2, 2, and 60Hz") lines(t, z1, col = 'red') lines(t, z2, col = 'blue') lines(t, z3, col = 'green') legend( "topleft", c("Input", "Pass: 0.11Hz", "Pass 15Hz", "Pass 1080Hz"), col = c(par("fg"), "red", "blue", "green"), lty = 1, cex = 0.6 ) # plot pwelch pwelch(x, fs = 2000, window = 4000, noverlap = 2000, plot = 1) pwelch(z1, fs = 2000, window = 4000, noverlap = 2000, plot = 2, col = "red") pwelch(z2, fs = 2000, window = 4000, noverlap = 2000, plot = 2, col = "blue") pwelch(z3, fs = 2000, window = 4000, noverlap = 2000, plot = 2, col = "green") #  Clean this demo  par(oldpar)
Provides five methods to baseline an array and calculate contrast.
baseline_array(x, along_dim, unit_dims = seq_along(dim(x))[along_dim], ...) ## S3 method for class 'array' baseline_array( x, along_dim, unit_dims = seq_along(dim(x))[along_dim], method = c("percentage", "sqrt_percentage", "decibel", "zscore", "sqrt_zscore", "subtract_mean"), baseline_indexpoints = NULL, baseline_subarray = NULL, ... )
baseline_array(x, along_dim, unit_dims = seq_along(dim(x))[along_dim], ...) ## S3 method for class 'array' baseline_array( x, along_dim, unit_dims = seq_along(dim(x))[along_dim], method = c("percentage", "sqrt_percentage", "decibel", "zscore", "sqrt_zscore", "subtract_mean"), baseline_indexpoints = NULL, baseline_subarray = NULL, ... )
x 
array (tensor) to calculate contrast 
along_dim 
integer range from 1 to the maximum dimension of 
unit_dims 
integer vector, baseline unit: see Details. 
... 
passed to other methods 
method 
character, baseline method options are:

baseline_indexpoints 
integer vector, which index points are counted
into baseline window? Each index ranges from 1 to 
baseline_subarray 
subarrays that should be used to calculate
baseline; default is 
Consider a scenario where we want to baseline a bunch of signals recorded
from different locations. For each location, we record n
sessions.
For each session, the signal is further decomposed into frequencytime
domain. In this case, we have the input x
in the following form:
$session x frequency x time x location$
Now we want to calibrate signals for each session, frequency and location using the first 100 time points as baseline points, then the code will be
$baseline_array(x, along_dim=3, baseline_window=1:100, unit_dims=c(1,2,4))$
along_dim=3
is dimension of time, in this case, it's the
third dimension of x
. baseline_indexpoints=1:100
, meaning
the first 100 time points are used to calculate baseline.
unit_dims
defines the unit signal. Its value c(1,2,4)
means the unit signal is per session (first dimension), per frequency
(second) and per location (fourth).
In some other cases, we might want to calculate baseline across frequencies
then the unit signal is $frequency x time$
, i.e. signals that share the
same session and location also share the same baseline. In this case,
we assign unit_dims=c(1,4)
.
There are five baseline methods. They fit for different types of data.
Denote $z$
is an unit signal, $z_0$
is its baseline slice. Then
these baseline methods are:
"percentage"
$\frac{z  \bar{z_{0}}}{\bar{z_{0}}} \times 100\%$
"sqrt_percentage"
$\frac{\sqrt{z}  \bar{\sqrt{z_{0}}}}{\bar{\sqrt{z_{0}}}} \times 100\%$
"decibel"
$10 \times ( \log_{10}(z)  \bar{\log_{10}(z_{0})} )$
"zscore"
$\frac{z\bar{z_{0}}}{sd(z_{0})}$
"sqrt_zscore"
$\frac{\sqrt{z}\bar{\sqrt{z_{0}}}}{sd(\sqrt{z_{0}})}$
Contrast array with the same dimension as x
.
# Set ncores = 2 to comply to CRAN policy. Please don't run this line ravetools_threads(n_threads = 2L) library(ravetools) set.seed(1) # Generate sample data dims = c(10,20,30,2) x = array(rnorm(prod(dims))^2, dims) # Set baseline window to be arbitrary 10 timepoints baseline_window = sample(30, 10) #  baseline percentage change  # Using base functions re1 < aperm(apply(x, c(1,2,4), function(y){ m < mean(y[baseline_window]) (y/m  1) * 100 }), c(2,3,1,4)) # Using ravetools re2 < baseline_array(x, 3, c(1,2,4), baseline_indexpoints = baseline_window, method = 'percentage') # Check different, should be very tiny (double precisions) range(re2  re1) # Check speed for large dataset, might take a while to profile ravetools_threads(n_threads = 1) dims < c(200,20,300,2) x < array(rnorm(prod(dims))^2, dims) # Set baseline window to be arbitrary 10 timepoints baseline_window < seq_len(100) f1 < function(){ aperm(apply(x, c(1,2,4), function(y){ m < mean(y[baseline_window]) (y/m  1) * 100 }), c(2,3,1,4)) } f2 < function(){ # equivalent as bl = x[,,baseline_window, ] # baseline_array(x, along_dim = 3, baseline_indexpoints = baseline_window, unit_dims = c(1,2,4), method = 'percentage') } range(f1()  f2()) microbenchmark::microbenchmark(f1(), f2(), times = 10L)
# Set ncores = 2 to comply to CRAN policy. Please don't run this line ravetools_threads(n_threads = 2L) library(ravetools) set.seed(1) # Generate sample data dims = c(10,20,30,2) x = array(rnorm(prod(dims))^2, dims) # Set baseline window to be arbitrary 10 timepoints baseline_window = sample(30, 10) #  baseline percentage change  # Using base functions re1 < aperm(apply(x, c(1,2,4), function(y){ m < mean(y[baseline_window]) (y/m  1) * 100 }), c(2,3,1,4)) # Using ravetools re2 < baseline_array(x, 3, c(1,2,4), baseline_indexpoints = baseline_window, method = 'percentage') # Check different, should be very tiny (double precisions) range(re2  re1) # Check speed for large dataset, might take a while to profile ravetools_threads(n_threads = 1) dims < c(200,20,300,2) x < array(rnorm(prod(dims))^2, dims) # Set baseline window to be arbitrary 10 timepoints baseline_window < seq_len(100) f1 < function(){ aperm(apply(x, c(1,2,4), function(y){ m < mean(y[baseline_window]) (y/m  1) * 100 }), c(2,3,1,4)) } f2 < function(){ # equivalent as bl = x[,,baseline_window, ] # baseline_array(x, along_dim = 3, baseline_indexpoints = baseline_window, unit_dims = c(1,2,4), method = 'percentage') } range(f1()  f2()) microbenchmark::microbenchmark(f1(), f2(), times = 10L)
Large filter order might not be optimal, but at lease this function
provides a feasible upper bound for the order such that the
filter has a stable AR
component.
butter_max_order( w, type = c("low", "high", "pass", "stop"), r = 10 * log10(2), tol = .Machine$double.eps )
butter_max_order( w, type = c("low", "high", "pass", "stop"), r = 10 * log10(2), tol = .Machine$double.eps )
w 
scaled frequency ranging from 0 to 1, where 1 is 'Nyquist' frequency 
type 
filter type 
r 
decibel attenuation at frequency 
tol 
tolerance of reciprocal condition number, default is

'Butterworth' filter in 'Arma' form.
# Find highest order (sharpest transition) of a bandpass filter sample_rate < 500 nyquist < sample_rate / 2 type < "pass" w < c(1, 50) / nyquist Rs < 6 # power attenuation at w # max order filter filter < butter_max_order(w, "pass", Rs) # 6 dB cutoff should be around 1 ~ 50 Hz diagnose_filter(filter$b, filter$a, fs = sample_rate)
# Find highest order (sharpest transition) of a bandpass filter sample_rate < 500 nyquist < sample_rate / 2 type < "pass" w < c(1, 50) / nyquist Rs < 6 # power attenuation at w # max order filter filter < butter_max_order(w, "pass", Rs) # 6 dB cutoff should be around 1 ~ 50 Hz diagnose_filter(filter$b, filter$a, fs = sample_rate)
Check 'Arma' filter
check_filter(b, a, w = NULL, r_expected = NULL, fs = NULL)
check_filter(b, a, w = NULL, r_expected = NULL, fs = NULL)
b 
moving average ( 
a 
autoregressive ( 
w 
normalized frequency, ranging from 0 to 1, where 1 is 'Nyquist' 
r_expected 
attenuation in decibel of each 
fs 
sample rate, used to infer the frequencies and formatting print message, not used in calculation; leave it blank by default 
A list of power estimation and the reciprocal condition number
of the AR
coefficients.
# create a butterworth filter with 3dB (halfpower) at [1, 5] Hz # and 60dB stopband attenuation at [0.5, 6] Hz sample_rate < 20 nyquist < sample_rate / 2 specs < buttord( Wp = c(1, 5) / nyquist, Ws = c(0.5, 6) / nyquist, Rp = 3, Rs = 60 ) filter < butter(specs) # filter quality is poor because the ARcoefficients # creates singular matrix with unstable inverse, # this will cause `filtfilt` to fail check_filter( b = filter$b, a = filter$a, # frequencies (normalized) where power is evaluated w = c(1, 5, 0.5, 6) / nyquist, # expected power r_expected = c(3, 3, 60, 60) )
# create a butterworth filter with 3dB (halfpower) at [1, 5] Hz # and 60dB stopband attenuation at [0.5, 6] Hz sample_rate < 20 nyquist < sample_rate / 2 specs < buttord( Wp = c(1, 5) / nyquist, Ws = c(0.5, 6) / nyquist, Rp = 3, Rs = 60 ) filter < butter(specs) # filter quality is poor because the ARcoefficients # creates singular matrix with unstable inverse, # this will cause `filtfilt` to fail check_filter( b = filter$b, a = filter$a, # frequencies (normalized) where power is evaluated w = c(1, 5, 0.5, 6) / nyquist, # expected power r_expected = c(3, 3, 60, 60) )
Collapse array
collapse(x, keep, ...) ## S3 method for class 'array' collapse( x, keep, average = TRUE, transform = c("asis", "10log10", "square", "sqrt"), ... )
collapse(x, keep, ...) ## S3 method for class 'array' collapse( x, keep, average = TRUE, transform = c("asis", "10log10", "square", "sqrt"), ... )
x 
A numeric multimode tensor (array), without 
keep 
Which dimension to keep 
... 
passed to other methods 
average 
collapse to sum or mean 
transform 
transform on the data before applying collapsing;
choices are 
a collapsed array with values to be mean or summation along collapsing dimensions
# Set ncores = 2 to comply to CRAN policy. Please don't run this line ravetools_threads(n_threads = 2L) # Example 1 x = matrix(1:16, 4) # Keep the first dimension and calculate sums along the rest collapse(x, keep = 1) rowMeans(x) # Should yield the same result # Example 2 x = array(1:120, dim = c(2,3,4,5)) result = collapse(x, keep = c(3,2)) compare = apply(x, c(3,2), mean) sum(abs(result  compare)) # The same, yield 0 or very small number (1e10) ravetools_threads(n_threads = 1) # Example 3 (performance) # Small data, no big difference x = array(rnorm(240), dim = c(4,5,6,2)) microbenchmark::microbenchmark( result = collapse(x, keep = c(3,2)), compare = apply(x, c(3,2), mean), times = 1L, check = function(v){ max(abs(range(do.call('', v)))) < 1e10 } ) # large data big difference x = array(rnorm(prod(300,200,105)), c(300,200,105,1)) microbenchmark::microbenchmark( result = collapse(x, keep = c(3,2)), compare = apply(x, c(3,2), mean), times = 1L , check = function(v){ max(abs(range(do.call('', v)))) < 1e10 })
# Set ncores = 2 to comply to CRAN policy. Please don't run this line ravetools_threads(n_threads = 2L) # Example 1 x = matrix(1:16, 4) # Keep the first dimension and calculate sums along the rest collapse(x, keep = 1) rowMeans(x) # Should yield the same result # Example 2 x = array(1:120, dim = c(2,3,4,5)) result = collapse(x, keep = c(3,2)) compare = apply(x, c(3,2), mean) sum(abs(result  compare)) # The same, yield 0 or very small number (1e10) ravetools_threads(n_threads = 1) # Example 3 (performance) # Small data, no big difference x = array(rnorm(240), dim = c(4,5,6,2)) microbenchmark::microbenchmark( result = collapse(x, keep = c(3,2)), compare = apply(x, c(3,2), mean), times = 1L, check = function(v){ max(abs(range(do.call('', v)))) < 1e10 } ) # large data big difference x = array(rnorm(prod(300,200,105)), c(300,200,105,1)) microbenchmark::microbenchmark( result = collapse(x, keep = c(3,2)), compare = apply(x, c(3,2), mean), times = 1L , check = function(v){ max(abs(range(do.call('', v)))) < 1e10 })
1D
, 2D
, 3D
data via FFT
Use the 'FastFourier' transform to compute the convolutions of two data
with zero padding. This function is mainly designed for image convolution.
For forward and backward convolution/filter, see filtfilt
.
convolve_signal(x, filter) convolve_image(x, filter) convolve_volume(x, filter)
convolve_signal(x, filter) convolve_image(x, filter) convolve_volume(x, filter)
x 
onedimensional signal vector, twodimensional image, or threedimensional volume; numeric or complex 
filter 
kernel with the same number of dimensions as 
This implementation uses 'FastFourier' transform to perform
1D
, 2D
, or 3D
convolution. Compared to implementations
using original mathematical definition of convolution, this approach is
much faster, especially for image and volume convolutions.
The input x
is zeropadded beyond edges. This is most common in image
or volume convolution, but less optimal for periodic onedimensional signals.
Please use other implementations if nonzero padding is needed.
The convolution results might be different to the ground truth by a precision
error, usually at 1e13
level, depending on the 'FFTW3'
library precision and implementation.
Convolution results with the same length and dimensions as x
.
If x
is complex, results will be complex, otherwise results will
be real numbers.
#  1D convolution  x < cumsum(rnorm(100)) filter < dnorm(2:2) # normalize filter < filter / sum(filter) smoothed < convolve_signal(x, filter) plot(x, pch = 20) lines(smoothed, col = 'red') #  2D convolution  x < array(0, c(100, 100)) x[ floor(runif(10, min = 1, max = 100)), floor(runif(10, min = 1, max = 100)) ] < 1 # smooth kernel < outer(dnorm(2:2), dnorm(2:2), FUN = "*") kernel < kernel / sum(kernel) y < convolve_image(x, kernel) oldpar < par(mfrow = c(1,2)) image(x, asp = 1, axes = FALSE, main = "Origin") image(y, asp = 1, axes = FALSE, main = "Smoothed") par(oldpar)
#  1D convolution  x < cumsum(rnorm(100)) filter < dnorm(2:2) # normalize filter < filter / sum(filter) smoothed < convolve_signal(x, filter) plot(x, pch = 20) lines(smoothed, col = 'red') #  2D convolution  x < array(0, c(100, 100)) x[ floor(runif(10, min = 1, max = 100)), floor(runif(10, min = 1, max = 100)) ] < 1 # smooth kernel < outer(dnorm(2:2), dnorm(2:2), FUN = "*") kernel < kernel / sum(kernel) y < convolve_image(x, kernel) oldpar < par(mfrow = c(1,2)) image(x, asp = 1, axes = FALSE, main = "Origin") image(y, asp = 1, axes = FALSE, main = "Smoothed") par(oldpar)
Decimate with 'FIR' or 'IIR' filter
decimate(x, q, n = if (ftype == "iir") 8 else 30, ftype = "fir")
decimate(x, q, n = if (ftype == "iir") 8 else 30, ftype = "fir")
x 
signal to be decimated 
q 
integer factor to downsample by 
n 
filter order used in the downsampling; default is 
ftype 
filter type, choices are 
This function is migrated from gsignal
package,
but with padding and indexing fixed. The results agree with 'Matlab'.
Decimated signal
x < 1:100 y < decimate(x, 2, ftype = "fir") y # compare with signal package z < gsignal::decimate(x, 2, ftype = "fir") # Compare decimated results plot(x, type = 'l') points(seq(1,100, 2), y, col = "green") points(seq(1,100, 2), z, col = "red")
x < 1:100 y < decimate(x, 2, ftype = "fir") y # compare with signal package z < gsignal::decimate(x, 2, ftype = "fir") # Compare decimated results plot(x, type = 'l') points(seq(1,100, 2), y, col = "green") points(seq(1,100, 2), z, col = "red")
Provides 'FIR' and 'IIR' filter options; default is 'FIR', see also
design_filter_fir
; for 'IIR' filters, see
design_filter_iir
.
design_filter( sample_rate, data = NULL, method = c("fir_kaiser", "firls", "fir_remez", "butter", "cheby1", "cheby2", "ellip"), high_pass_freq = NA, high_pass_trans_freq = NA, low_pass_freq = NA, low_pass_trans_freq = NA, passband_ripple = 0.1, stopband_attenuation = 40, filter_order = NA, ..., data_size = length(data) )
design_filter( sample_rate, data = NULL, method = c("fir_kaiser", "firls", "fir_remez", "butter", "cheby1", "cheby2", "ellip"), high_pass_freq = NA, high_pass_trans_freq = NA, low_pass_freq = NA, low_pass_trans_freq = NA, passband_ripple = 0.1, stopband_attenuation = 40, filter_order = NA, ..., data_size = length(data) )
sample_rate 
data sample rate 
data 
data to be filtered, can be optional ( 
method 
filter method, options are 
high_pass_freq , low_pass_freq

highpass or lowpass frequency,
see 
high_pass_trans_freq , low_pass_trans_freq

transition bandwidths,
see 
passband_ripple 
allowable passband ripple in decibel; default is

stopband_attenuation 
minimum stopband attenuation (in decibel) at
transition frequency; default is 
filter_order 
suggested filter order; 'RAVE' may or may not adopt this suggestion depending on the data and numerical feasibility 
... 
passed to filter generator functions 
data_size 
used by 'FIR' filter design to determine maximum order,
ignored in 'IIR' filters; automatically derived from 
If data
is specified and nonempty, this function returns
filtered data via forward and backward filtfilt
; if data
is
NULL
, then returns the generator function.
sample_rate < 200 t < seq(0, 10, by = 1 / sample_rate) x < sin(t * 4 * pi) + sin(t * 20 * pi) + 2 * sin(t * 120 * pi) + rnorm(length(t), sd = 0.4) #  Using FIR  # Lowpass filter y1 < design_filter( data = x, sample_rate = sample_rate, low_pass_freq = 3, low_pass_trans_freq = 0.5 ) # Bandpass cheby1 filter 812 Hz with custom transition y2 < design_filter( data = x, method = "cheby1", sample_rate = sample_rate, low_pass_freq = 12, low_pass_trans_freq = .25, high_pass_freq = 8, high_pass_trans_freq = .25 ) y3 < design_filter( data = x, sample_rate = sample_rate, low_pass_freq = 80, high_pass_freq = 30 ) oldpar < par(mfrow = c(2, 1), mar = c(3.1, 2.1, 3.1, 0.1)) plot(t, x, type = 'l', xlab = "Time", ylab = "", main = "Mixture of 2, 10, and 60Hz", xlim = c(0,1)) # lines(t, y, col = 'red') lines(t, y3, col = 'green') lines(t, y2, col = 'blue') lines(t, y1, col = 'red') legend( "topleft", c("Input", "Low: 3Hz", "Pass 812Hz", "Pass 3080Hz"), col = c(par("fg"), "red", "blue", "green"), lty = 1, cex = 0.6 ) # plot pwelch pwelch(x, fs = sample_rate, window = sample_rate * 2, noverlap = sample_rate, plot = 1, ylim = c(100, 10)) pwelch(y1, fs = sample_rate, window = sample_rate * 2, noverlap = sample_rate, plot = 2, col = "red") pwelch(y2, fs = sample_rate, window = sample_rate * 2, noverlap = sample_rate, plot = 2, col = "blue") pwelch(y3, fs = sample_rate, window = sample_rate * 2, noverlap = sample_rate, plot = 2, col = "green") #  Clean this demo  par(oldpar)
sample_rate < 200 t < seq(0, 10, by = 1 / sample_rate) x < sin(t * 4 * pi) + sin(t * 20 * pi) + 2 * sin(t * 120 * pi) + rnorm(length(t), sd = 0.4) #  Using FIR  # Lowpass filter y1 < design_filter( data = x, sample_rate = sample_rate, low_pass_freq = 3, low_pass_trans_freq = 0.5 ) # Bandpass cheby1 filter 812 Hz with custom transition y2 < design_filter( data = x, method = "cheby1", sample_rate = sample_rate, low_pass_freq = 12, low_pass_trans_freq = .25, high_pass_freq = 8, high_pass_trans_freq = .25 ) y3 < design_filter( data = x, sample_rate = sample_rate, low_pass_freq = 80, high_pass_freq = 30 ) oldpar < par(mfrow = c(2, 1), mar = c(3.1, 2.1, 3.1, 0.1)) plot(t, x, type = 'l', xlab = "Time", ylab = "", main = "Mixture of 2, 10, and 60Hz", xlim = c(0,1)) # lines(t, y, col = 'red') lines(t, y3, col = 'green') lines(t, y2, col = 'blue') lines(t, y1, col = 'red') legend( "topleft", c("Input", "Low: 3Hz", "Pass 812Hz", "Pass 3080Hz"), col = c(par("fg"), "red", "blue", "green"), lty = 1, cex = 0.6 ) # plot pwelch pwelch(x, fs = sample_rate, window = sample_rate * 2, noverlap = sample_rate, plot = 1, ylim = c(100, 10)) pwelch(y1, fs = sample_rate, window = sample_rate * 2, noverlap = sample_rate, plot = 2, col = "red") pwelch(y2, fs = sample_rate, window = sample_rate * 2, noverlap = sample_rate, plot = 2, col = "blue") pwelch(y3, fs = sample_rate, window = sample_rate * 2, noverlap = sample_rate, plot = 2, col = "green") #  Clean this demo  par(oldpar)
'FIR'
filter using
firls
Design 'FIR'
filter using firls
design_filter_fir( sample_rate, filter_order = NA, data_size = NA, high_pass_freq = NA, high_pass_trans_freq = NA, low_pass_freq = NA, low_pass_trans_freq = NA, stopband_attenuation = 40, scale = TRUE, method = c("kaiser", "firls", "remez") )
design_filter_fir( sample_rate, filter_order = NA, data_size = NA, high_pass_freq = NA, high_pass_trans_freq = NA, low_pass_freq = NA, low_pass_trans_freq = NA, stopband_attenuation = 40, scale = TRUE, method = c("kaiser", "firls", "remez") )
sample_rate 
sampling frequency 
filter_order 
filter order, leave 
data_size 
minimum length of data to apply the filter, used to
decide the maximum filter order. For 'FIR' filter, data length must be
greater than 
high_pass_freq 
highpass frequency; default is 
high_pass_trans_freq 
highpass frequency bandwidth; default
is automatically inferred from data size.
Frequency 
low_pass_freq 
lowpass frequency; default is 
low_pass_trans_freq 
lowpass frequency bandwidth; default
is automatically inferred from data size.
Frequency 
stopband_attenuation 
allowable power attenuation (in decibel) at
transition frequency; default is 
scale 
whether to scale the filter for unity gain 
method 
method to generate 'FIR' filter, default is using

Filter type is determined from high_pass_freq
and
low_pass_freq
. Highpass frequency is ignored if high_pass_freq
is NA
, hence the filter is lowpass filter. When
low_pass_freq
is NA
, then
the filter is highpass filter. When both high_pass_freq
and
low_pass_freq
are valid (positive, less than 'Nyquist'), then
the filter is a bandpass filter if bandpass is less than lowpass
frequency, otherwise the filter is bandstop.
Although the peak amplitudes are set at 1 by low_pass_freq
and
high_pass_freq
, the transition from peak amplitude to zero require
a transition, which is tricky but also important to set.
When 'FIR' filters have too steep transition boundaries, the filter tends to
have ripples in peak amplitude, introducing artifacts to the final signals.
When the filter is too flat, components from unwanted frequencies may also
get aliased into the filtered signals. Ideally, the transition bandwidth
cannot be too steep nor too flat. In this function, users may control
the transition frequency bandwidths via low_pass_trans_freq
and
high_pass_trans_freq
. The power at the end of transition is defined
by stopband_attenuation
, with default value of 40
(i.e.
40 dB, this number is automatically negated during the calculation).
By design, a lowpass 5 Hz filter with 1 Hz transition bandwidth results in
around 40 dB power at 6 Hz.
'FIR' filter in 'Arma' form.
#  Basic  sample_rate < 500 data_size < 1000 # lowpass at 5 Hz, with auto transition bandwidth # from kaiser's method, with default stopband attenuation = 40 dB filter < design_filter_fir( low_pass_freq = 5, sample_rate = sample_rate, data_size = data_size ) # Passband ripple is around 0.08 dB # stopband attenuation is around 40 dB print(filter) diagnose_filter( filter$b, filter$a, fs = sample_rate, n = data_size, cutoffs = c(3, 6, 40), vlines = 5 ) #  Advanced  sample_rate < 500 data_size < 1000 # Rejecting 38 Hz, with transition bandwidth 0.5 Hz at both ends # Using leastsquare (firls) to generate FIR filter # Suggesting the filter order n=160 filter < design_filter_fir( low_pass_freq = 3, low_pass_trans_freq = 0.5, high_pass_freq = 8, high_pass_trans_freq = 0.5, filter_order = 160, sample_rate = sample_rate, data_size = data_size, method = "firls" ) # print(filter) diagnose_filter( filter$b, filter$a, fs = sample_rate, n = data_size, cutoffs = c(1, 40), vlines = c(3, 8) )
#  Basic  sample_rate < 500 data_size < 1000 # lowpass at 5 Hz, with auto transition bandwidth # from kaiser's method, with default stopband attenuation = 40 dB filter < design_filter_fir( low_pass_freq = 5, sample_rate = sample_rate, data_size = data_size ) # Passband ripple is around 0.08 dB # stopband attenuation is around 40 dB print(filter) diagnose_filter( filter$b, filter$a, fs = sample_rate, n = data_size, cutoffs = c(3, 6, 40), vlines = 5 ) #  Advanced  sample_rate < 500 data_size < 1000 # Rejecting 38 Hz, with transition bandwidth 0.5 Hz at both ends # Using leastsquare (firls) to generate FIR filter # Suggesting the filter order n=160 filter < design_filter_fir( low_pass_freq = 3, low_pass_trans_freq = 0.5, high_pass_freq = 8, high_pass_trans_freq = 0.5, filter_order = 160, sample_rate = sample_rate, data_size = data_size, method = "firls" ) # print(filter) diagnose_filter( filter$b, filter$a, fs = sample_rate, n = data_size, cutoffs = c(1, 40), vlines = c(3, 8) )
Design an 'IIR' filter
design_filter_iir( method = c("butter", "cheby1", "cheby2", "ellip"), sample_rate, filter_order = NA, high_pass_freq = NA, high_pass_trans_freq = NA, low_pass_freq = NA, low_pass_trans_freq = NA, passband_ripple = 0.1, stopband_attenuation = 40 )
design_filter_iir( method = c("butter", "cheby1", "cheby2", "ellip"), sample_rate, filter_order = NA, high_pass_freq = NA, high_pass_trans_freq = NA, low_pass_freq = NA, low_pass_trans_freq = NA, passband_ripple = 0.1, stopband_attenuation = 40 )
method 
filter method name, choices are 
sample_rate 
sampling frequency 
filter_order 
suggested filter order. Notice filters with higher orders
may become numerically unstable, hence this number is only a suggested
number. If the filter is unstable, this function will choose a lower order;
leave this input 
high_pass_freq 
highpass frequency; default is 
high_pass_trans_freq 
highpass frequency bandwidth; default is automatically inferred from filter type. 
low_pass_freq 
lowpass frequency; default is 
low_pass_trans_freq 
lowpass frequency bandwidth; default is automatically inferred from filter type. 
passband_ripple 
allowable passband ripple in decibel; default is

stopband_attenuation 
minimum stopband attenuation (in decibel) at
transition frequency; default is 
A filter in 'Arma' form.
sample_rate < 500 my_diagnose < function( filter, vlines = c(8, 12), cutoffs = c(3, 6)) { diagnose_filter( b = filter$b, a = filter$a, fs = sample_rate, vlines = vlines, cutoffs = cutoffs ) } #  Default using butterworth to generate 812 bandpass filter  # Butterworth filter with cutoff frequency # 7 ~ 13 (default transition bandwidth is 1Hz) at 3 dB filter < design_filter_iir( method = "butter", low_pass_freq = 12, high_pass_freq = 8, sample_rate = 500 ) filter my_diagnose(filter) ## explicit bandwidths and attenuation (sharper transition) # Butterworth filter with cutoff frequency # passband ripple is 0.5 dB (812 Hz) # stopband attenuation is 40 dB (518 Hz) filter < design_filter_iir( method = "butter", low_pass_freq = 12, low_pass_trans_freq = 6, high_pass_freq = 8, high_pass_trans_freq = 3, sample_rate = 500, passband_ripple = 0.5, stopband_attenuation = 40 ) filter my_diagnose(filter) #  cheby1  filter < design_filter_iir( method = "cheby1", low_pass_freq = 12, high_pass_freq = 8, sample_rate = 500 ) my_diagnose(filter) #  cheby2  filter < design_filter_iir( method = "cheby2", low_pass_freq = 12, high_pass_freq = 8, sample_rate = 500 ) my_diagnose(filter) #  ellip  filter < design_filter_iir( method = "ellip", low_pass_freq = 12, high_pass_freq = 8, sample_rate = 500 ) my_diagnose(filter)
sample_rate < 500 my_diagnose < function( filter, vlines = c(8, 12), cutoffs = c(3, 6)) { diagnose_filter( b = filter$b, a = filter$a, fs = sample_rate, vlines = vlines, cutoffs = cutoffs ) } #  Default using butterworth to generate 812 bandpass filter  # Butterworth filter with cutoff frequency # 7 ~ 13 (default transition bandwidth is 1Hz) at 3 dB filter < design_filter_iir( method = "butter", low_pass_freq = 12, high_pass_freq = 8, sample_rate = 500 ) filter my_diagnose(filter) ## explicit bandwidths and attenuation (sharper transition) # Butterworth filter with cutoff frequency # passband ripple is 0.5 dB (812 Hz) # stopband attenuation is 40 dB (518 Hz) filter < design_filter_iir( method = "butter", low_pass_freq = 12, low_pass_trans_freq = 6, high_pass_freq = 8, high_pass_trans_freq = 3, sample_rate = 500, passband_ripple = 0.5, stopband_attenuation = 40 ) filter my_diagnose(filter) #  cheby1  filter < design_filter_iir( method = "cheby1", low_pass_freq = 12, high_pass_freq = 8, sample_rate = 500 ) my_diagnose(filter) #  cheby2  filter < design_filter_iir( method = "cheby2", low_pass_freq = 12, high_pass_freq = 8, sample_rate = 500 ) my_diagnose(filter) #  ellip  filter < design_filter_iir( method = "ellip", low_pass_freq = 12, high_pass_freq = 8, sample_rate = 500 ) my_diagnose(filter)
'Detrending' is often used before the signal power calculation.
detrend(x, trend = c("constant", "linear"), break_points = NULL)
detrend(x, trend = c("constant", "linear"), break_points = NULL)
x 
numerical or complex, a vector or a matrix 
trend 
the trend of the signal; choices are 
break_points 
integer vector, or 
The signals with trend removed in matrix form; the number of columns is the number of signals, and number of rows is length of the signals
x < rnorm(100, mean = 1) + c( seq(0, 5, length.out = 50), seq(5, 3, length.out = 50)) plot(x) plot(detrend(x, 'constant')) plot(detrend(x, 'linear')) plot(detrend(x, 'linear', 50))
x < rnorm(100, mean = 1) + c( seq(0, 5, length.out = 50), seq(5, 3, length.out = 50)) plot(x) plot(detrend(x, 'constant')) plot(detrend(x, 'linear')) plot(detrend(x, 'linear', 50))
The diagnostic plots include 'Welch Periodogram'
(pwelch
) and histogram (hist
)
diagnose_channel( s1, s2 = NULL, sc = NULL, srate, name = "", try_compress = TRUE, max_freq = 300, window = ceiling(srate * 2), noverlap = window/2, std = 3, which = NULL, main = "Channel Inspection", col = c("black", "red"), cex = 1.2, cex.lab = 1, lwd = 0.5, plim = NULL, nclass = 100, start_time = 0, boundary = NULL, mar = c(3.1, 4.1, 2.1, 0.8) * (0.25 + cex * 0.75) + 0.1, mgp = cex * c(2, 0.5, 0), xaxs = "i", yaxs = "i", xline = 1.66 * cex, yline = 2.66 * cex, tck = 0.005 * (3 + cex), ... )
diagnose_channel( s1, s2 = NULL, sc = NULL, srate, name = "", try_compress = TRUE, max_freq = 300, window = ceiling(srate * 2), noverlap = window/2, std = 3, which = NULL, main = "Channel Inspection", col = c("black", "red"), cex = 1.2, cex.lab = 1, lwd = 0.5, plim = NULL, nclass = 100, start_time = 0, boundary = NULL, mar = c(3.1, 4.1, 2.1, 0.8) * (0.25 + cex * 0.75) + 0.1, mgp = cex * c(2, 0.5, 0), xaxs = "i", yaxs = "i", xline = 1.66 * cex, yline = 2.66 * cex, tck = 0.005 * (3 + cex), ... )
s1 
the main signal to draw 
s2 
the comparing signal to draw; usually 
sc 
decimated 
srate 
sampling rate 
name 
name of 
try_compress 
whether try to compress (decimate) 
max_freq 
the maximum frequency to display in 'Welch Periodograms' 
window , noverlap

see 
std 
the standard deviation of the channel signals used to determine

which 

main 
the title of the signal plot 
col 
colors of 
cex , lwd , mar , cex.lab , mgp , xaxs , yaxs , tck , ...

graphical parameters; see

plim 
the yaxis limit to draw in 'Welch Periodograms' 
nclass 
number of classes to show in histogram
( 
start_time 
the starting time of channel (will only be used to draw signals) 
boundary 
a red boundary to show in channel plot; default is
to be automatically determined by 
xline , yline

distance of axis labels towards ticks 
A list of boundary and yaxis limit used to draw the channel
library(ravetools) # Generate 20 second data at 2000 Hz time < seq(0, 20, by = 1 / 2000) signal < sin( 120 * pi * time) + sin(time * 20*pi) + exp(time^2) * cos(time * 10*pi) + rnorm(length(time)) signal2 < notch_filter(signal, 2000) diagnose_channel(signal, signal2, srate = 2000, name = c("Raw", "Filtered"), cex = 1)
library(ravetools) # Generate 20 second data at 2000 Hz time < seq(0, 20, by = 1 / 2000) signal < sin( 120 * pi * time) + sin(time * 20*pi) + exp(time^2) * cos(time * 10*pi) + rnorm(length(time)) signal2 < notch_filter(signal, 2000) diagnose_channel(signal, signal2, srate = 2000, name = c("Raw", "Filtered"), cex = 1)
Generate frequency response plot with sampledata simulation
diagnose_filter( b, a, fs, n = 512, whole = FALSE, sample = stats::rnorm(n, mean = sample_signal(n), sd = 0.2), vlines = NULL, xlim = "auto", cutoffs = c(3, 6, 12) )
diagnose_filter( b, a, fs, n = 512, whole = FALSE, sample = stats::rnorm(n, mean = sample_signal(n), sd = 0.2), vlines = NULL, xlim = "auto", cutoffs = c(3, 6, 12) )
b 
the movingaverage coefficients of an 
a 
the autoregressive coefficients of an 
fs 
sampling frequency in 
n 
number of points at which to evaluate the frequency response;
default is 
whole 
whether to evaluate beyond 
sample 
sample signal of length 
vlines 
additional vertical lines (frequencies) to plot 
xlim 
frequency limit of frequency response plot; default is

cutoffs 
cutoff decibel powers to draw on the frequency plot, also used
to calculate the frequency limit when 
Nothing
library(ravetools) # sample rate srate < 500 # signal length npts < 1000 # bandpass bpass < c(1, 50) # Nyquist fn < srate / 2 w < bpass / fn #  FIR filter  order < 160 # FIR1 is MA filter, a = 1 filter < fir1(order, w, "pass") diagnose_filter( b = filter$b, a = filter$a, n = npts, fs = srate, vlines = bpass ) #  Butter filter  filter < butter(3, w, "pass") diagnose_filter( b = filter$b, a = filter$a, n = npts, fs = srate, vlines = bpass )
library(ravetools) # sample rate srate < 500 # signal length npts < 1000 # bandpass bpass < c(1, 50) # Nyquist fn < srate / 2 w < bpass / fn #  FIR filter  order < 160 # FIR1 is MA filter, a = 1 filter < fir1(order, w, "pass") diagnose_filter( b = filter$b, a = filter$a, n = npts, fs = srate, vlines = bpass ) #  Butter filter  filter < butter(3, w, "pass") diagnose_filter( b = filter$b, a = filter$a, n = npts, fs = srate, vlines = bpass )
Calculate surface distances of graph or mesh using 'Dijkstra' method.
dijkstras_surface_distance( positions, faces, start_node, face_index_start = NA, max_search_distance = NA, ... ) surface_path(x, target_node)
dijkstras_surface_distance( positions, faces, start_node, face_index_start = NA, max_search_distance = NA, ... ) surface_path(x, target_node)
positions 
numeric matrix with no 
faces 
integer matrix with each row containing indices of nodes. For
graphs, 
start_node 
integer, row index of 
face_index_start 
integer, the start of the nodes in 
max_search_distance 
numeric, maximum distance to iterate;
default is 
... 
reserved for backward compatibility 
x 
distance calculation results returned by

target_node 
the target node number to reach (from the starting node);

dijkstras_surface_distance
returns a list distance
table with the meta configurations. surface_path
returns a data frame
of the node ID (from start_node
to target_node
) and cumulative
distance along the shortest path.
#  Toy example  # Position is 2D, total 6 points positions < matrix(runif(6 * 2), ncol = 2) # edges defines connected nodes edges < matrix(ncol = 2, byrow = TRUE, data = c( 1,2, 2,3, 1,3, 2,4, 3,4, 2,5, 4,5, 2,5, 4,6, 5,6 )) # calculate distances ret < dijkstras_surface_distance( start_node = 1, positions = positions, faces = edges, face_index_start = 1 ) # get shortest path from the first node to the last path < surface_path(ret, target_node = 6) # plot the results from_node < path$path[nrow(path)] to_node < path$path[1] plot(positions, pch = 16, axes = FALSE, xlab = "X", ylab = "Y", main = "Dijkstra's shortest path") segments( x0 = positions[edges[,1],1], y0 = positions[edges[,1],2], x1 = positions[edges[,2],1], y1 = positions[edges[,2],2] ) points(positions[path$path,], col = "steelblue", pch = 16) arrows( x0 = positions[from_node,1], y0 = positions[from_node,2], x1 = positions[to_node,1], y1 = positions[to_node,2], col = "steelblue", lwd = 2, length = 0.1, lty = 2 ) points(positions[1,,drop=FALSE], pch = 16, col = "orangered") points(positions[6,,drop=FALSE], pch = 16, col = "purple3") #  Example with mesh  ## Not run: # Please install the downstream package `threeBrain` # and call library(threeBrain) # the following code set up the files read.fs.surface < internal_rave_function( "read.fs.surface", "threeBrain") default_template_directory < internal_rave_function( "default_template_directory", "threeBrain") surface_path < file.path(default_template_directory(), "N27", "surf", "lh.pial") if(!file.exists(surface_path)) { internal_rave_function( "download_N27", "threeBrain")() } # Example starts from here > # Load the mesh mesh < read.fs.surface(surface_path) # Calculate the path with maximum radius 100 ret < dijkstras_surface_distance( start_node = 1, positions = mesh$vertices, faces = mesh$faces, max_search_distance = 100, verbose = TRUE ) # get shortest path from the first node to node 43144 path < surface_path(ret, target_node = 43144) # plot from_nodes < path$path[nrow(path)] to_nodes < path$path[1] # calculate colors pal < colorRampPalette( colors = c("red", "orange", "orange3", "purple3", "purple4") )(1001) col < pal[ceiling( path$distance / max(path$distance, na.rm = TRUE) * 1000 ) + 1] oldpar < par(mfrow = c(2, 2), mar = c(0, 0, 0, 0)) for(xdim in c(1, 2, 3)) { if( xdim < 3 ) { ydim < xdim + 1 } else { ydim < 3 xdim < 1 } plot( mesh$vertices[, xdim], mesh$vertices[, ydim], pch = ".", col = "#BEBEBE33", axes = FALSE, xlab = "P  A", ylab = "S  I", asp = 1 ) segments( x0 = mesh$vertices[from_nodes, xdim], y0 = mesh$vertices[from_nodes, ydim], x1 = mesh$vertices[to_nodes, xdim], y1 = mesh$vertices[to_nodes, ydim], col = col ) } # plot distance map distances < ret$paths$distance col < pal[ceiling(distances / max(distances, na.rm = TRUE) * 1000) + 1] selection < !is.na(distances) plot( mesh$vertices[, 2], mesh$vertices[, 3], pch = ".", col = "#BEBEBE33", axes = FALSE, xlab = "P  A", ylab = "S  I", asp = 1 ) points( mesh$vertices[selection, c(2, 3)], col = col[selection], pch = "." ) # reset graphic state par(oldpar) ## End(Not run)
#  Toy example  # Position is 2D, total 6 points positions < matrix(runif(6 * 2), ncol = 2) # edges defines connected nodes edges < matrix(ncol = 2, byrow = TRUE, data = c( 1,2, 2,3, 1,3, 2,4, 3,4, 2,5, 4,5, 2,5, 4,6, 5,6 )) # calculate distances ret < dijkstras_surface_distance( start_node = 1, positions = positions, faces = edges, face_index_start = 1 ) # get shortest path from the first node to the last path < surface_path(ret, target_node = 6) # plot the results from_node < path$path[nrow(path)] to_node < path$path[1] plot(positions, pch = 16, axes = FALSE, xlab = "X", ylab = "Y", main = "Dijkstra's shortest path") segments( x0 = positions[edges[,1],1], y0 = positions[edges[,1],2], x1 = positions[edges[,2],1], y1 = positions[edges[,2],2] ) points(positions[path$path,], col = "steelblue", pch = 16) arrows( x0 = positions[from_node,1], y0 = positions[from_node,2], x1 = positions[to_node,1], y1 = positions[to_node,2], col = "steelblue", lwd = 2, length = 0.1, lty = 2 ) points(positions[1,,drop=FALSE], pch = 16, col = "orangered") points(positions[6,,drop=FALSE], pch = 16, col = "purple3") #  Example with mesh  ## Not run: # Please install the downstream package `threeBrain` # and call library(threeBrain) # the following code set up the files read.fs.surface < internal_rave_function( "read.fs.surface", "threeBrain") default_template_directory < internal_rave_function( "default_template_directory", "threeBrain") surface_path < file.path(default_template_directory(), "N27", "surf", "lh.pial") if(!file.exists(surface_path)) { internal_rave_function( "download_N27", "threeBrain")() } # Example starts from here > # Load the mesh mesh < read.fs.surface(surface_path) # Calculate the path with maximum radius 100 ret < dijkstras_surface_distance( start_node = 1, positions = mesh$vertices, faces = mesh$faces, max_search_distance = 100, verbose = TRUE ) # get shortest path from the first node to node 43144 path < surface_path(ret, target_node = 43144) # plot from_nodes < path$path[nrow(path)] to_nodes < path$path[1] # calculate colors pal < colorRampPalette( colors = c("red", "orange", "orange3", "purple3", "purple4") )(1001) col < pal[ceiling( path$distance / max(path$distance, na.rm = TRUE) * 1000 ) + 1] oldpar < par(mfrow = c(2, 2), mar = c(0, 0, 0, 0)) for(xdim in c(1, 2, 3)) { if( xdim < 3 ) { ydim < xdim + 1 } else { ydim < 3 xdim < 1 } plot( mesh$vertices[, xdim], mesh$vertices[, ydim], pch = ".", col = "#BEBEBE33", axes = FALSE, xlab = "P  A", ylab = "S  I", asp = 1 ) segments( x0 = mesh$vertices[from_nodes, xdim], y0 = mesh$vertices[from_nodes, ydim], x1 = mesh$vertices[to_nodes, xdim], y1 = mesh$vertices[to_nodes, ydim], col = col ) } # plot distance map distances < ret$paths$distance col < pal[ceiling(distances / max(distances, na.rm = TRUE) * 1000) + 1] selection < !is.na(distances) plot( mesh$vertices[, 2], mesh$vertices[, 3], pch = ".", col = "#BEBEBE33", axes = FALSE, xlab = "P  A", ylab = "S  I", asp = 1 ) points( mesh$vertices[selection, c(2, 3)], col = col[selection], pch = "." ) # reset graphic state par(oldpar) ## End(Not run)
Speed up covariance calculation for large matrices. The
default behavior is the same as cov
('pearson'
,
no NA
handling).
fast_cov(x, y = NULL, col_x = NULL, col_y = NULL, df = NA)
fast_cov(x, y = NULL, col_x = NULL, col_y = NULL, df = NA)
x 
a numeric vector, matrix or data frame; a matrix is highly recommended to maximize the performance 
y 
NULL (default) or a vector, matrix or data frame with compatible
dimensions to x; the default is equivalent to 
col_x 
integers indicating the subset indices (columns) of 
col_y 
integers indicating the subset indices (columns) of 
df 
a scalar indicating the degrees of freedom; default is

A covariance matrix of x
and y
. Note that there is no
NA
handling. Any missing values will lead to NA
in the
resulting covariance matrices.
# Set ncores = 2 to comply to CRAN policy. Please don't run this line ravetools_threads(n_threads = 2L) x < matrix(rnorm(400), nrow = 100) # Call `cov(x)` to compare fast_cov(x) # Calculate covariance of subsets fast_cov(x, col_x = 1, col_y = 1:2) # Speed comparison, better to use multiple cores (4, 8, or more) # to show the differences. ravetools_threads(n_threads = 1) x < matrix(rnorm(100000), nrow = 1000) microbenchmark::microbenchmark( fast_cov = { fast_cov(x, col_x = 1:50, col_y = 51:100) }, cov = { cov(x[,1:50], x[,51:100]) }, unit = 'ms', times = 10 )
# Set ncores = 2 to comply to CRAN policy. Please don't run this line ravetools_threads(n_threads = 2L) x < matrix(rnorm(400), nrow = 100) # Call `cov(x)` to compare fast_cov(x) # Calculate covariance of subsets fast_cov(x, col_x = 1, col_y = 1:2) # Speed comparison, better to use multiple cores (4, 8, or more) # to show the differences. ravetools_threads(n_threads = 1) x < matrix(rnorm(100000), nrow = 1000) microbenchmark::microbenchmark( fast_cov = { fast_cov(x, col_x = 1:50, col_y = 51:100) }, cov = { cov(x[,1:50], x[,51:100]) }, unit = 'ms', times = 10 )
Compute quantiles
fast_quantile(x, prob = 0.5, na.rm = FALSE, ...) fast_median(x, na.rm = FALSE, ...) fast_mvquantile(x, prob = 0.5, na.rm = FALSE, ...) fast_mvmedian(x, na.rm = FALSE, ...)
fast_quantile(x, prob = 0.5, na.rm = FALSE, ...) fast_median(x, na.rm = FALSE, ...) fast_mvquantile(x, prob = 0.5, na.rm = FALSE, ...) fast_mvmedian(x, na.rm = FALSE, ...)
x 
numericalvalue vector for 
prob 
a probability with value from 0 to 1 
na.rm 
logical; if true, any 
... 
reserved for future use 
fast_quantile
and fast_median
calculate univariate
quantiles (singlevalue return); fast_mvquantile
and fast_mvmedian
calculate multivariate quantiles (for each column, result lengths equal to
the number of columns).
fast_quantile(runif(1000), 0.1) fast_median(1:100) x < matrix(rnorm(100), ncol = 2) fast_mvquantile(x, 0.2) fast_mvmedian(x) # Compare speed for vectors (usually 30% faster) x < rnorm(10000) microbenchmark::microbenchmark( fast_median = fast_median(x), base_median = median(x), # bioc_median = Biobase::rowMedians(matrix(x, nrow = 1)), times = 100, unit = "milliseconds" ) # Multivariate cases # (5~7x faster than base R) # (3~5x faster than Biobase rowMedians) x < matrix(rnorm(100000), ncol = 20) microbenchmark::microbenchmark( fast_median = fast_mvmedian(x), base_median = apply(x, 2, median), # bioc_median = Biobase::rowMedians(t(x)), times = 10, unit = "milliseconds" )
fast_quantile(runif(1000), 0.1) fast_median(1:100) x < matrix(rnorm(100), ncol = 2) fast_mvquantile(x, 0.2) fast_mvmedian(x) # Compare speed for vectors (usually 30% faster) x < rnorm(10000) microbenchmark::microbenchmark( fast_median = fast_median(x), base_median = median(x), # bioc_median = Biobase::rowMedians(matrix(x, nrow = 1)), times = 100, unit = "milliseconds" ) # Multivariate cases # (5~7x faster than base R) # (3~5x faster than Biobase rowMedians) x < matrix(rnorm(100000), ncol = 20) microbenchmark::microbenchmark( fast_median = fast_mvmedian(x), base_median = apply(x, 2, median), # bioc_median = Biobase::rowMedians(t(x)), times = 10, unit = "milliseconds" )
Create a cube volume (256
'voxels' on each margin), fill
in the 'voxels' that are inside of the surface.
fill_surface( surface, inflate = 0, IJK2RAS = NULL, preview = FALSE, preview_frame = 128 )
fill_surface( surface, inflate = 0, IJK2RAS = NULL, preview = FALSE, preview_frame = 128 )
surface 
a surface mesh, can be mesh objects from 
inflate 
amount of 'voxels' to inflate on the final result; must be
a nonnegative integer. A zero 
IJK2RAS 
volume 'IJK' (zeroindexed coordinate index) to

preview 
whether to preview the results; default is false 
preview_frame 
integer from 1 to 256 the depth frame used to generate preview. 
This function creates a volume (256 on each margin) and fill in the volume from a surface mesh. The surface vertex points will be embedded into the volume first. These points may not be connected together, hence for each 'voxel', a cube patch will be applied to grow the volume. Then, the volume will be bucketfilled from a corner, forming a negated mask of "outsideofsurface" area. The inverted bucketfilled volume is then shrunk so the mask boundary tightly fits the surface
A list containing the filled volume and parameters used to generate the volume
Zhengjia Wang
# takes > 5s to run example # Generate a sphere surface < vcg_sphere() surface$vb[1:3, ] < surface$vb[1:3, ] * 50 fill_surface(surface, preview = TRUE)
# takes > 5s to run example # Generate a sphere surface < vcg_sphere() surface$vb[1:3, ] < surface$vb[1:3, ] * 50 fill_surface(surface, preview = TRUE)
The function is written from the scratch. The result has been
compared against the 'Matlab' filter
function with onedimensional
real inputs. Other situations such as matrix b
or multidimensional
x
are not implemented. For double filters (forwardbackward),
see filtfilt
.
filter_signal(b, a, x, z)
filter_signal(b, a, x, z)
b 
onedimensional real numerical vector, the movingaverage
coefficients of an 
a 
the autoregressive (recursive) coefficients of an 
x 
numerical vector input (real value) 
z 
initial condition, must have length of 
A list of two vectors: the first vector is the filtered signal;
the second vector is the final state of z
t < seq(0, 1, by = 0.01) x < sin(2 * pi * t * 2.3) bf < gsignal::butter(2, c(0.15, 0.3)) res < filter_signal(bf$b, bf$a, x) y < res[[1]] z < res[[2]] ## Matlab (2022a) equivalent: # t = [0:0.01:1]; # x = sin(2 * pi * t * 2.3); # [b,a] = butter(2,[.15,.3]); # [y,z] = filter(b, a, x)
t < seq(0, 1, by = 0.01) x < sin(2 * pi * t * 2.3) bf < gsignal::butter(2, c(0.15, 0.3)) res < filter_signal(bf$b, bf$a, x) y < res[[1]] z < res[[2]] ## Matlab (2022a) equivalent: # t = [0:0.01:1]; # x = sin(2 * pi * t * 2.3); # [b,a] = butter(2,[.15,.3]); # [y,z] = filter(b, a, x)
Filter window functions
hanning(n) hamming(n) blackman(n) blackmannuttall(n) blackmanharris(n) flattopwin(n) bohmanwin(n)
hanning(n) hamming(n) blackman(n) blackmannuttall(n) blackmanharris(n) flattopwin(n) bohmanwin(n)
n 
number of timepoints in window 
A numeric vector of window with length n
hanning(10) hamming(11) blackmanharris(21)
hanning(10) hamming(11) blackmanharris(21)
The result has been tested against 'Matlab' filtfilt
function. Currently this function only supports one filter at a time.
filtfilt(b, a, x)
filtfilt(b, a, x)
b 
onedimensional real numerical vector, the movingaverage
coefficients of an 
a 
the autoregressive (recursive) coefficients of an 
x 
numerical vector input (real value) 
The filtered signal, normally the same length as the input signal
x
.
t < seq(0, 1, by = 0.01) x < sin(2 * pi * t * 2.3) bf < gsignal::butter(2, c(0.15, 0.3)) res < filtfilt(bf$b, bf$a, x) ## Matlab (2022a) equivalent: # t = [0:0.01:1]; # x = sin(2 * pi * t * 2.3); # [b,a] = butter(2,[.15,.3]); # res = filtfilt(b, a, x)
t < seq(0, 1, by = 0.01) x < sin(2 * pi * t * 2.3) bf < gsignal::butter(2, c(0.15, 0.3)) res < filtfilt(bf$b, bf$a, x) ## Matlab (2022a) equivalent: # t = [0:0.01:1]; # x = sin(2 * pi * t * 2.3); # [b,a] = butter(2,[.15,.3]); # res = filtfilt(b, a, x)
FIR
filter designGenerate a fir1
filter that is checked against Matlab
fir1
function.
fir1( n, w, type = c("low", "high", "stop", "pass", "DC0", "DC1"), window = hamming, scale = TRUE, hilbert = FALSE )
fir1( n, w, type = c("low", "high", "stop", "pass", "DC0", "DC1"), window = hamming, scale = TRUE, hilbert = FALSE )
n 
filter order 
w 
band edges, nondecreasing vector in the range 0 to 1, where 1 is
the 
type 
type of the filter, one of 
window 
smoothing window function or a numerical vector. The filter is
the same shape as the smoothing window. When 
scale 
whether to scale the filter; default is true 
hilbert 
whether to use 'Hilbert' transformer; default is false 
The FIR
filter coefficients with class 'Arma'
.
The moving average coefficient is a vector of length n+1
.
FIR
filter designProduce a linear phase filter from the weighted mean squared such that error in the specified bands is minimized.
firls(N, freq, A, W = NULL, ftype = "")
firls(N, freq, A, W = NULL, ftype = "")
N 
filter order, must be even (if odd, then will be increased by one) 
freq 
vector of frequency points in the range from 0 to 1, where 1
corresponds to the 
A 
vector of the same length as 
W 
weighting function that contains one value for each band that
weights the mean squared error in that band. 
ftype 
transformer type; default is 
The FIR
filter coefficients with class 'Arma'
.
The moving average coefficient is a vector of length n+1
.
Compute the zplane frequency response of an ARMA
model.
freqz2(b, a = 1, fs = 2 * pi, n = 512, whole = FALSE, ...)
freqz2(b, a = 1, fs = 2 * pi, n = 512, whole = FALSE, ...)
b 
the movingaverage coefficients of an 
a 
the autoregressive coefficients of an 
fs 
sampling frequency in 
n 
number of points at which to evaluate the frequency response;
default is 
whole 
whether to evaluate beyond 
... 
ignored 
A list of frequencies and corresponding responses in complex vector
Grow volume mask
grow_volume(volume, x, y = x, z = x, threshold = 0.5)
grow_volume(volume, x, y = x, z = x, threshold = 0.5)
volume 
volume mask array, must be 3dimensional array 
x , y , z

size of grow along each direction 
threshold 
threshold after convolution 
A binary volume mask
oldpar < par(mfrow = c(2,3), mar = c(0.1,0.1,3.1,0.1)) mask < array(0, c(21,21,21)) mask[11,11,11] < 1 image(mask[11,,], asp = 1, main = "Original mask", axes = FALSE) image(grow_volume(mask, 2)[11,,], asp = 1, main = "Dilated (size=2) mask", axes = FALSE) image(grow_volume(mask, 5)[11,,], asp = 1, main = "Dilated (size=5) mask", axes = FALSE) mask[11, sample(11,2), sample(11,2)] < 1 image(mask[11,,], asp = 1, main = "Original mask", axes = FALSE) image(grow_volume(mask, 2)[11,,], asp = 1, main = "Dilated (size=2) mask", axes = FALSE) image(grow_volume(mask, 5)[11,,], asp = 1, main = "Dilated (size=5) mask", axes = FALSE) par(oldpar)
oldpar < par(mfrow = c(2,3), mar = c(0.1,0.1,3.1,0.1)) mask < array(0, c(21,21,21)) mask[11,11,11] < 1 image(mask[11,,], asp = 1, main = "Original mask", axes = FALSE) image(grow_volume(mask, 2)[11,,], asp = 1, main = "Dilated (size=2) mask", axes = FALSE) image(grow_volume(mask, 5)[11,,], asp = 1, main = "Dilated (size=5) mask", axes = FALSE) mask[11, sample(11,2), sample(11,2)] < 1 image(mask[11,,], asp = 1, main = "Original mask", axes = FALSE) image(grow_volume(mask, 2)[11,,], asp = 1, main = "Dilated (size=2) mask", axes = FALSE) image(grow_volume(mask, 5)[11,,], asp = 1, main = "Dilated (size=5) mask", axes = FALSE) par(oldpar)
Internal function used for examples relative to 'RAVE' project and should not be used directly.
internal_rave_function(name, pkg, inherit = TRUE, on_missing = NULL)
internal_rave_function(name, pkg, inherit = TRUE, on_missing = NULL)
name 
function or variable name 
pkg 
'RAVE' package name 
inherit 
passed to 
on_missing 
default value to return of no function is found 
Function object if found, otherwise on_missing
.
Find and interpolate stimulation signals
interpolate_stimulation( x, sample_rate, duration = 40/sample_rate, ord = 4L, nknots = 100, nsd = 1, nstim = NULL, regularization = 0.5 )
interpolate_stimulation( x, sample_rate, duration = 40/sample_rate, ord = 4L, nknots = 100, nsd = 1, nstim = NULL, regularization = 0.5 )
x 
numerical vector representing a analog signal 
sample_rate 
sampling frequency 
duration 
time in second: duration of interpolation 
ord 
spline order, default is 4 
nknots 
a rough number of knots to use, default is 100 
nsd 
number of standard deviation to detect stimulation signals, default is 1 
nstim 
number of stimulation pulses, default is to autodetect 
regularization 
regularization parameter in case of inverting singular matrices, default is 0.5 
Interpolated signal with an attribute of which sample points are interpolated
x0 < rnorm(1000) / 5 + sin(1:1000 / 300) # Simulates pulase signals x < x0 x[400:410] < 100 x[420:430] < 100 fitted < interpolate_stimulation(x, 100, duration = 0.3, nknots = 10, nsd = 2) oldpar < par(mfrow = c(2, 1)) plot(fitted, type = 'l', col = 'blue', lwd = 2) lines(x, col = 'red') lines(x0, col = 'black') legend("topleft", c("Interpolated", "Observed", "Underlying"), lty = 1, col = c("blue", "red", "black")) pwelch(x0, 100, 200, 100, plot = 1, col = 'black', ylim = c(50, 50)) pwelch(x, 100, 200, 100, plot = 2, col = 'red') pwelch(fitted, 100, 200, 100, plot = 2, col = 'blue') par(oldpar)
x0 < rnorm(1000) / 5 + sin(1:1000 / 300) # Simulates pulase signals x < x0 x[400:410] < 100 x[420:430] < 100 fitted < interpolate_stimulation(x, 100, duration = 0.3, nknots = 10, nsd = 2) oldpar < par(mfrow = c(2, 1)) plot(fitted, type = 'l', col = 'blue', lwd = 2) lines(x, col = 'red') lines(x0, col = 'black') legend("topleft", c("Interpolated", "Observed", "Underlying"), lty = 1, col = c("blue", "red", "black")) pwelch(x0, 100, 200, 100, plot = 1, col = 'black', ylim = c(50, 50)) pwelch(x, 100, 200, 100, plot = 2, col = 'red') pwelch(fitted, 100, 200, 100, plot = 2, col = 'blue') par(oldpar)
Left 'Hippocampus' of 'N27Collin' brain
left_hippocampus_mask
left_hippocampus_mask
A threemode integer mask array with values of 1
('Hippocampus')
and 0
(other brain tissues)
'Matlab' heatmap plot palette
matlab_palette()
matlab_palette()
vector of 64 colors
This function is softdeprecated. Please use
vcg_mesh_volume
, vcg_uniform_remesh
, and
vcg_smooth_explicit
or vcg_smooth_implicit
.
mesh_from_volume( volume, output_format = c("rgl", "freesurfer"), IJK2RAS = NULL, threshold = 0, verbose = TRUE, remesh = TRUE, remesh_voxel_size = 1, remesh_multisample = TRUE, remesh_automerge = TRUE, smooth = FALSE, smooth_lambda = 10, smooth_delta = 20, smooth_method = "surfPreserveLaplace" )
mesh_from_volume( volume, output_format = c("rgl", "freesurfer"), IJK2RAS = NULL, threshold = 0, verbose = TRUE, remesh = TRUE, remesh_voxel_size = 1, remesh_multisample = TRUE, remesh_automerge = TRUE, smooth = FALSE, smooth_lambda = 10, smooth_delta = 20, smooth_method = "surfPreserveLaplace" )
volume 
3dimensional volume array 
output_format 
resulting data format, choices are 
IJK2RAS 
volume 'IJK' (zeroindexed coordinate index) to

threshold 
threshold used to create volume mask; the surface will be created to fit the mask boundaries 
verbose 
whether to verbose the progress 
remesh 
whether to resample the mesh using 
remesh_voxel_size , remesh_multisample , remesh_automerge

see
arguments in 
smooth 
whether to smooth the mesh via 
smooth_lambda , smooth_delta , smooth_method

A 'mesh3d'
surface if output_format
is 'rgl', or
'fs.surface'
surface otherwise.
volume < array(0, dim = c(8,8,8)) volume[4:5, 4:5, 4:5] < 1 graphics::image(x = volume[4,,]) # you can use rgl::wire3d(mesh) to visualize the mesh mesh < mesh_from_volume(volume, verbose = FALSE)
volume < array(0, dim = c(8,8,8)) volume[4:5, 4:5, 4:5] < 1 graphics::image(x = volume[4,,]) # you can use rgl::wire3d(mesh) to visualize the mesh mesh < mesh_from_volume(volume, verbose = FALSE)
Compute 'multitaper' spectral densities of timeseries data
multitaper_config( data_length, fs, frequency_range = NULL, time_bandwidth = 5, num_tapers = NULL, window_params = c(5, 1), nfft = NA, detrend_opt = "linear" ) multitaper( data, fs, frequency_range = NULL, time_bandwidth = 5, num_tapers = NULL, window_params = c(5, 1), nfft = NA, detrend_opt = "linear" )
multitaper_config( data_length, fs, frequency_range = NULL, time_bandwidth = 5, num_tapers = NULL, window_params = c(5, 1), nfft = NA, detrend_opt = "linear" ) multitaper( data, fs, frequency_range = NULL, time_bandwidth = 5, num_tapers = NULL, window_params = c(5, 1), nfft = NA, detrend_opt = "linear" )
data_length 
length of data 
fs 
sampling frequency in 'Hz' 
frequency_range 
frequency range to look at; length of two 
time_bandwidth 
a number indicating timehalf bandwidth product; i.e.
the window duration times the half bandwidth of main lobe; default is

num_tapers 
number of 'DPSS' tapers to use; default is 
window_params 
vector of two numbers; the first number is the
window size in seconds; the second number if the step size; default is

nfft 
'NFFT' size, positive; see 'Details' 
detrend_opt 
how you want to remove the trend from data window; options
are 
data 
numerical vector, signal traces 
The original source code comes from 'Prerau' Lab (see 'Github'
repository 'multitaper_toolbox'
under user 'preraulab'
).
The results tend to agree with their 'Python' implementation with precision
on the order of at 1E7
with standard deviation at most 1E5
.
The original copy was licensed under a Creative Commons Attribution
'NC''SA' 4.0 International License
(https://creativecommons.org/licenses/byncsa/4.0/).
This package ('ravetools'
) redistributes the multitaper
function under minor modifications on nfft
. In the original copy
there is no parameter to control the exact numbers of nfft
, and
the nfft
is always the power of 2. While choosing
nfft
to be the power of 2 is always recommended, the modified code
allows other choices.
multitaper_config
returns a list of configuration parameters
for the filters; multitaper
also returns the time, frequency and
corresponding spectral power.
# Takes long to run time < seq(0, 3, by = 0.001) x < sin(time * 20*pi) + exp(time^2) * cos(time * 10*pi) res < multitaper( x, 1000, frequency_range = c(0,15), time_bandwidth=1.5, window_params=c(2,0.01) ) image( x = res$time, y = res$frequency, z = 10 * log10(res$spec), xlab = "Time (s)", ylab = 'Frequency (Hz)', col = matlab_palette() )
# Takes long to run time < seq(0, 3, by = 0.001) x < sin(time * 20*pi) + exp(time^2) * cos(time * 10*pi) res < multitaper( x, 1000, frequency_range = c(0,15), time_bandwidth=1.5, window_params=c(2,0.01) ) image( x = res$time, y = res$frequency, z = 10 * log10(res$spec), xlab = "Time (s)", ylab = 'Frequency (Hz)', col = matlab_palette() )
Matrix4
instance for 'Affine'
transformCreate a Matrix4
instance for 'Affine'
transform
new_matrix4() as_matrix4(m)
new_matrix4() as_matrix4(m)
m 
a matrix or a vector to be converted to the 
A Matrix4
instance
Quaternion
instance to store '3D' rotationCreate instances that mimic the 'three.js'
syntax.
new_quaternion(x = 0, y = 0, z = 0, w = 1) as_quaternion(q)
new_quaternion(x = 0, y = 0, z = 0, w = 1) as_quaternion(q)
x , y , z , w

numeric of length one 
q 
R object to be converted to 
A Quaternion
instance
Vector3
instance to store '3D' pointsCreate instances that mimic the 'three.js'
syntax.
new_vector3(x = 0, y = 0, z = 0) as_vector3(v)
new_vector3(x = 0, y = 0, z = 0) as_vector3(v)
x , y , z

numeric, must have the same length, 
v 
R object to be converted to 
A Vector3
instance
vec3 < new_vector3( x = 1:9, y = 9:1, z = rep(c(1,2,3), 3) ) vec3[] # transform m < new_matrix4() # rotation xy plane by 30 degrees m$make_rotation_z(pi / 6) vec3$apply_matrix4(m) vec3[] as_vector3(c(1,2,3))
vec3 < new_vector3( x = 1:9, y = 9:1, z = rep(c(1,2,3), 3) ) vec3[] # transform m < new_matrix4() # rotation xy plane by 30 degrees m$make_rotation_z(pi / 6) vec3$apply_matrix4(m) vec3[] as_vector3(c(1,2,3))
Apply 'Notch' filter
notch_filter( s, sample_rate, lb = c(59, 118, 178), ub = c(61, 122, 182), domain = 1 )
notch_filter( s, sample_rate, lb = c(59, 118, 178), ub = c(61, 122, 182), domain = 1 )
s 
numerical vector if 
sample_rate 
sample rate 
lb 
filter lower bound of the frequencies to remove 
ub 
filter upper bound of the frequencies to remove;
shares the same length as 
domain 

Mainly used to remove electrical line frequencies
at 60, 120, and 180 Hz
.
filtered signal in time domain (real numerical vector)
time < seq(0, 3, 0.005) s < sin(120 * pi * time) + rnorm(length(time)) # Welch periodogram shows a peak at 60Hz pwelch(s, 200, plot = 1, log = "y") # notch filter to remove 60Hz s1 < notch_filter(s, 200, lb = 59, ub = 61) pwelch(s1, 200, plot = 2, log = "y", col = "red")
time < seq(0, 3, 0.005) s < sin(120 * pi * time) + rnorm(length(time)) # Welch periodogram shows a peak at 60Hz pwelch(s, 200, plot = 1, log = "y") # notch filter to remove 60Hz s1 < notch_filter(s, 200, lb = 59, ub = 61) pwelch(s1, 200, plot = 2, log = "y", col = "red")
Set or get thread options
detect_threads() ravetools_threads(n_threads = "auto", stack_size = "auto")
detect_threads() ravetools_threads(n_threads = "auto", stack_size = "auto")
n_threads 
number of threads to set 
stack_size 
Stack size (in bytes) to use for worker threads. The
default used for 
detect_threads
returns an integer of default threads that
is determined by the number of CPU
cores; ravetools_threads
returns nothing.
detect_threads() ravetools_threads(n_threads = 2)
detect_threads() ravetools_threads(n_threads = 2)
Plot one or more signal traces in the same figure
plot_signals( signals, sample_rate = 1, col = graphics::par("fg"), space = 0.995, space_mode = c("quantile", "absolute"), start_time = 0, duration = NULL, compress = TRUE, channel_names = NULL, time_shift = 0, xlab = "Time (s)", ylab = "Electrode", lwd = 0.5, new_plot = TRUE, xlim = NULL, cex = 1, cex.lab = 1, mar = c(3.1, 2.1, 2.1, 0.8) * (0.25 + cex * 0.75) + 0.1, mgp = cex * c(2, 0.5, 0), xaxs = "r", yaxs = "i", xline = 1.5 * cex, yline = 1 * cex, tck = 0.005 * (3 + cex), ... )
plot_signals( signals, sample_rate = 1, col = graphics::par("fg"), space = 0.995, space_mode = c("quantile", "absolute"), start_time = 0, duration = NULL, compress = TRUE, channel_names = NULL, time_shift = 0, xlab = "Time (s)", ylab = "Electrode", lwd = 0.5, new_plot = TRUE, xlim = NULL, cex = 1, cex.lab = 1, mar = c(3.1, 2.1, 2.1, 0.8) * (0.25 + cex * 0.75) + 0.1, mgp = cex * c(2, 0.5, 0), xaxs = "r", yaxs = "i", xline = 1.5 * cex, yline = 1 * cex, tck = 0.005 * (3 + cex), ... )
signals 
numerical matrix with each row to be a signal trace and each column contains the signal values at a time point 
sample_rate 
sampling frequency 
col 
signal color, can be vector of one or more 
space 
vertical spacing among the traces; for values greater than 1,
the spacing is absolute; default is 
space_mode 
mode of spacing, only used when 
start_time 
the time to start drawing relative to the first column 
duration 
duration of the signal to draw 
compress 
whether to compress signals if the data is too large 
channel_names 

time_shift 
the actual start time of the signal. Unlike

xlab , ylab , lwd , xlim , cex , cex.lab , mar , mgp , xaxs , yaxs , tck , ...


new_plot 
whether to draw a new plot; default is true 
xline , yline

the gap between axis and label 
n < 1000 base_signal < c(rep(0, n/2), sin(seq(0,10,length.out = n/2))) * 10 signals < rbind(rnorm(n) + base_signal, rbinom(n, 10, 0.3) + base_signal, rt(n, 5) + base_signal) plot_signals(signals, sample_rate = 100) plot_signals(signals, sample_rate = 100, start_time = 5) plot_signals(signals, sample_rate = 100, start_time = 5, time_shift = 100)
n < 1000 base_signal < c(rep(0, n/2), sin(seq(0,10,length.out = n/2))) * 10 signals < rbind(rnorm(n) + base_signal, rbinom(n, 10, 0.3) + base_signal, rt(n, 5) + base_signal) plot_signals(signals, sample_rate = 100) plot_signals(signals, sample_rate = 100, start_time = 5) plot_signals(signals, sample_rate = 100, start_time = 5, time_shift = 100)
pwelch
is for single signal trace only; mv_pwelch
is for multiple traces. Currently mv_pwelch
is experimental and
should not be called directly.
pwelch( x, fs, window = 64, noverlap = window/2, nfft = "auto", window_family = hamming, col = "black", xlim = NULL, ylim = NULL, main = "Welch periodogram", plot = 0, log = c("xy", "", "x", "y"), ... ) ## S3 method for class ''ravetoolspwelch'' print(x, ...) ## S3 method for class ''ravetoolspwelch'' plot( x, log = c("xy", "x", "y", ""), se = FALSE, xticks, type = "l", add = FALSE, col = graphics::par("fg"), col.se = "orange", alpha.se = 0.5, lty = 1, lwd = 1, cex = 1, las = 1, main = "Welch periodogram", xlab, ylab, xlim = NULL, ylim = NULL, xaxs = "i", yaxs = "i", xline = 1.2 * cex, yline = 2 * cex, mar = c(2.6, 3.8, 2.1, 0.6) * (0.5 + cex/2), mgp = cex * c(2, 0.5, 0), tck = 0.02 * cex, grid = TRUE, ... ) mv_pwelch( x, margin, fs, window = 64, noverlap = window/2, nfft = "auto", window_family = hamming )
pwelch( x, fs, window = 64, noverlap = window/2, nfft = "auto", window_family = hamming, col = "black", xlim = NULL, ylim = NULL, main = "Welch periodogram", plot = 0, log = c("xy", "", "x", "y"), ... ) ## S3 method for class ''ravetoolspwelch'' print(x, ...) ## S3 method for class ''ravetoolspwelch'' plot( x, log = c("xy", "x", "y", ""), se = FALSE, xticks, type = "l", add = FALSE, col = graphics::par("fg"), col.se = "orange", alpha.se = 0.5, lty = 1, lwd = 1, cex = 1, las = 1, main = "Welch periodogram", xlab, ylab, xlim = NULL, ylim = NULL, xaxs = "i", yaxs = "i", xline = 1.2 * cex, yline = 2 * cex, mar = c(2.6, 3.8, 2.1, 0.6) * (0.5 + cex/2), mgp = cex * c(2, 0.5, 0), tck = 0.02 * cex, grid = TRUE, ... ) mv_pwelch( x, margin, fs, window = 64, noverlap = window/2, nfft = "auto", window_family = hamming )
x 
numerical vector or a rowmajor vector, signals.
If 
fs 
sample rate, average number of time points per second 
window 
window length in time points, default size is 
noverlap 
overlap between two adjacent windows, measured in time
points; default is half of the 
nfft 
number of points in window function; default is automatically determined from input data and window, scaled up to the nearest power of 2 
window_family 
function generator for generating filter windows,
default is 
col , xlim , ylim , main , type , cex , las , xlab , ylab , lty , lwd , xaxs , yaxs , mar , mgp , tck

parameters passed to 
plot 
integer, whether to plot the result or not; choices are 
log 
indicates which axis should be 
... 
will be passed to 
se 
logical or a positive number indicating whether to plot standard error of mean; default is false. If provided with a number, then a multiple of standard error will be drawn. This option is only available when power is in logscale (decibel unit) 
xticks 
ticks to show on frequency axis 
add 
logical, whether the plot should be added to existing canvas 
col.se , alpha.se

controls the color and opacity of the standard error 
xline , yline

controls how close the axis labels to the corresponding axes 
grid 
whether to draw rectangular grid lines to the plot; only
respected when 
margin 
the margin in which 
A list with class 'ravetoolspwelch'
that contains the
following items:
freq
frequencies used to calculate the 'periodogram'
spec
resulting spectral power for each frequency
window
window function (in numerical vector) used
noverlap
number of overlapping timepoints between two adjacent windows
nfft
number of basis functions
fs
sample rate
x_len
input signal length
method
a character string 'Welch'
x < rnorm(1000) pwel < pwelch(x, 100) pwel plot(pwel, log = "xy")
x < rnorm(1000) pwel < pwelch(x, 100) pwel plot(pwel, log = "xy")
Convert raw vectors to R vectors
raw_to_uint8(x) raw_to_uint16(x) raw_to_uint32(x) raw_to_int8(x) raw_to_int16(x) raw_to_int32(x) raw_to_int64(x) raw_to_float(x) raw_to_string(x)
raw_to_uint8(x) raw_to_uint16(x) raw_to_uint32(x) raw_to_int8(x) raw_to_int16(x) raw_to_int32(x) raw_to_int64(x) raw_to_float(x) raw_to_string(x)
x 
raw vector of bytes 
For numeric conversions, the function names are straightforward.
For example,
raw_to_uintN
converts raw vectors to unsigned integers, and
raw_to_intN
converts raw vectors to signed integers. The number
'N'
stands for the number of bits used to store the integer.
For example raw_to_uint8
uses 8 bits (1 byte) to store an integer,
hence the value range is 0255
.
The input data length must be multiple of the element size represented by
the underlying data. For example uint16
integer uses 16 bites, and
one raw number uses 8 bits, hence two raw vectors can form one unsigned
integer16. That is, raw_to_uint16
requires the length of input
to be multiple of two. An easy calculation is: the length of x
times
8, must be divided by 'N'
(see last paragraph for definition).
The returned data uses the closest available R native data type that can
fully represent the data. For example, R does not have single float
type, hence raw_to_float
returns double
type, which can
represent all possible values in float
. For raw_to_uint32
,
the potential value range is 0  (2^321)
. This exceeds the limit of
R integer type (2^31)  (2^311)
. Therefore, the returned values
will be real (double float) data type.
There is no native data type that can store integer64 data in R, package
bit64
provides integer64
type, which will be used by
raw_to_int64
. Currently there is no solution to convert raw to
unsigned integer64 type.
raw_to_string
converts raw to character string. This function respects
null
character, hence is slightly different than the native
rawToChar
, which translates raw bytebybyte. If each
raw byte represents a valid character, then the above two functions returns
the same result. However, when the characters represented by raw bytes are
invalid, raw_to_string
will stop parsing and returns only the valid
characters, while rawToChar
will still try to parse, and
most likely to result in errors.
Please see Examples for comparisons.
Numeric vectors, except for raw_to_string
, which returns
a string.
# 0x00, 0x7f, 0x80, 0xFF x < as.raw(c(0, 127, 128, 255)) raw_to_uint8(x) # The first bit becomes the integer sign # 128 > 128, 255 > 1 raw_to_int8(x) ## Comments based on little endian system # 0x7f00 (32512), 0xFF80 (65408 unsigned, or 128 signed) raw_to_uint16(x) raw_to_int16(x) # 0xFF807F00 (4286611200 unsigned, 8356096 signed) raw_to_uint32(x) raw_to_int32(x) #  String  # ASCII case: all valid x < charToRaw("This is an ASCII string") raw_to_string(x) rawToChar(x) x < c(charToRaw("This is the end."), as.raw(0), charToRaw("*** is invalid")) # rawToChar will raise error raw_to_string(x) #  Integer64  # Runs on little endian system x < as.raw(c(0x80, 0x00, 0x7f, 0x80, 0xFF, 0x50, 0x7f, 0x00)) # Calculate bitstring, which concaternates the followings # 10000000 (0x80), 00000000 (0x00), 01111111 (0x7f), 10000000 (0x80), # 11111111 (0xFF), 01010000 (0x50), 01111111 (0x7f), 00000000 (0x00) if(.Platform$endian == "little") { bitstring < paste0( "00000000011111110101000011111111", "10000000011111110000000010000000" ) } else { bitstring < paste0( "00000001000000001111111000000001", "11111111000010101111111000000000" ) } # This is expected value bit64::as.integer64(structure( bitstring, class = "bitstring" )) # This is actual value raw_to_int64(x)
# 0x00, 0x7f, 0x80, 0xFF x < as.raw(c(0, 127, 128, 255)) raw_to_uint8(x) # The first bit becomes the integer sign # 128 > 128, 255 > 1 raw_to_int8(x) ## Comments based on little endian system # 0x7f00 (32512), 0xFF80 (65408 unsigned, or 128 signed) raw_to_uint16(x) raw_to_int16(x) # 0xFF807F00 (4286611200 unsigned, 8356096 signed) raw_to_uint32(x) raw_to_int32(x) #  String  # ASCII case: all valid x < charToRaw("This is an ASCII string") raw_to_string(x) rawToChar(x) x < c(charToRaw("This is the end."), as.raw(0), charToRaw("*** is invalid")) # rawToChar will raise error raw_to_string(x) #  Integer64  # Runs on little endian system x < as.raw(c(0x80, 0x00, 0x7f, 0x80, 0xFF, 0x50, 0x7f, 0x00)) # Calculate bitstring, which concaternates the followings # 10000000 (0x80), 00000000 (0x00), 01111111 (0x7f), 10000000 (0x80), # 11111111 (0xFF), 01010000 (0x50), 01111111 (0x7f), 00000000 (0x00) if(.Platform$endian == "little") { bitstring < paste0( "00000000011111110101000011111111", "10000000011111110000000010000000" ) } else { bitstring < paste0( "00000001000000001111111000000001", "11111111000010101111111000000000" ) } # This is expected value bit64::as.integer64(structure( bitstring, class = "bitstring" )) # This is actual value raw_to_int64(x)
Test whether the filter is numerically stable for filtfilt
.
rcond_filter_ar(a)
rcond_filter_ar(a)
a 
autoregression coefficient, numerical vector; the first element must not be zero 
Reciprocal condition number of matrix z1
, used in
filtfilt
. If the number is less than
.Machine$double.eps
, then filtfilt
will fail.
# Butterworth filter with lowpass at 0.1 Hz (order = 4) filter < butter(4, 0.1, "low") # TRUE rcond_filter_ar(filter$a) > .Machine$double.eps diagnose_filter(filter$b, filter$a, 500) # Bad filter (order is too high) filter < butter(50, 0.1, "low") rcond_filter_ar(filter$a) > .Machine$double.eps # filtfilt needs to inverse a singular matrix diagnose_filter(filter$b, filter$a, 500)
# Butterworth filter with lowpass at 0.1 Hz (order = 4) filter < butter(4, 0.1, "low") # TRUE rcond_filter_ar(filter$a) > .Machine$double.eps diagnose_filter(filter$b, filter$a, 500) # Bad filter (order is too high) filter < butter(50, 0.1, "low") rcond_filter_ar(filter$a) > .Machine$double.eps # filtfilt needs to inverse a singular matrix diagnose_filter(filter$b, filter$a, 500)
'NiftyReg'
Registers 'CT' to 'MRI', or 'MRI' to another 'MRI'
register_volume( source, target, method = c("rigid", "affine", "nonlinear"), interpolation = c("cubic", "trilinear", "nearest"), threads = detect_threads(), symmetric = TRUE, verbose = TRUE, ... )
register_volume( source, target, method = c("rigid", "affine", "nonlinear"), interpolation = c("cubic", "trilinear", "nearest"), threads = detect_threads(), symmetric = TRUE, verbose = TRUE, ... )
source 
source imaging data, or a 
target 
target imaging data to align to; for example, 'MRI' 
method 
method of transformation, choices are 
interpolation 
how volumes should be interpolated, choices are

threads , symmetric , verbose , ...

see 
See niftyreg
source < system.file("extdata", "epi_t2.nii.gz", package="RNiftyReg") target < system.file("extdata", "flash_t1.nii.gz", package="RNiftyReg") aligned < register_volume(source, target, verbose = FALSE) source_img < aligned$source[[1]] target_img < aligned$target aligned_img < aligned$image oldpar < par(mfrow = c(2, 2), mar = c(0.1, 0.1, 3.1, 0.1)) pal < grDevices::grey.colors(256, alpha = 1) image(source_img[,,30], asp = 1, axes = FALSE, col = pal, main = "Source image") image(target_img[,,64], asp = 1, axes = FALSE, col = pal, main = "Target image") image(aligned_img[,,64], asp = 1, axes = FALSE, col = pal, main = "Aligned image") # bucket fill and calculate differences aligned_img[is.nan(aligned_img)  aligned_img <= 1] < 1 target_img[is.nan(target_img)  aligned_img <= 1] < 1 diff < abs(aligned_img / target_img  1) image(diff[,,64], asp = 1, axes = FALSE, col = pal, main = "Percentage Difference") par(oldpar)
source < system.file("extdata", "epi_t2.nii.gz", package="RNiftyReg") target < system.file("extdata", "flash_t1.nii.gz", package="RNiftyReg") aligned < register_volume(source, target, verbose = FALSE) source_img < aligned$source[[1]] target_img < aligned$target aligned_img < aligned$image oldpar < par(mfrow = c(2, 2), mar = c(0.1, 0.1, 3.1, 0.1)) pal < grDevices::grey.colors(256, alpha = 1) image(source_img[,,30], asp = 1, axes = FALSE, col = pal, main = "Source image") image(target_img[,,64], asp = 1, axes = FALSE, col = pal, main = "Target image") image(aligned_img[,,64], asp = 1, axes = FALSE, col = pal, main = "Aligned image") # bucket fill and calculate differences aligned_img[is.nan(aligned_img)  aligned_img <= 1] < 1 target_img[is.nan(target_img)  aligned_img <= 1] < 1 diff < abs(aligned_img / target_img  1) image(diff[,,64], asp = 1, axes = FALSE, col = pal, main = "Percentage Difference") par(oldpar)
'rgl'
without requiring 'x11'
Internally used for example showcases. Please install package 'rgl'
manually to use these functions.
rgl_call(FUN, ...) rgl_view(expr, quoted = FALSE, env = parent.frame()) rgl_plot_normals(x, length = 1, lwd = 1, col = 1, ...)
rgl_call(FUN, ...) rgl_view(expr, quoted = FALSE, env = parent.frame()) rgl_plot_normals(x, length = 1, lwd = 1, col = 1, ...)
FUN 

... 
passed to 
expr 
expression within which 
quoted 
whether 
env 
environment in which 
x 
triangular 
length , lwd , col

normal vector length, size, and color 
# Make sure the example does not run when compiling # or check the package if(FALSE) { volume < array(0, dim = c(8,8,8)) volume[4:5, 4:5, 4:5] < 1 mesh < mesh_from_volume(volume, verbose = FALSE) rgl_view({ rgl_call("shade3d", mesh, col = 3) rgl_plot_normals(mesh) }) }
# Make sure the example does not run when compiling # or check the package if(FALSE) { volume < array(0, dim = c(8,8,8)) volume[4:5, 4:5, 4:5] < 1 mesh < mesh_from_volume(volume, verbose = FALSE) rgl_view({ rgl_call("shade3d", mesh, col = 3) rgl_plot_normals(mesh) }) }
Rearrange arrays in parallel
shift_array(x, along_margin, unit_margin, shift_amount)
shift_array(x, along_margin, unit_margin, shift_amount)
x 
array, must have at least matrix 
along_margin 
which index is to be shifted 
unit_margin 
which dimension decides 
shift_amount 
shift amount along 
A simple usecase for this function is to think of a matrix where each row is a signal and columns stand for time. The objective is to align (timelock) each signal according to certain events. For each signal, we want to shift the time points by certain amount.
In this case, the shift amount is defined by shift_amount
, whose
length equals to number of signals. along_margin=2
as we want to shift
time points (column, the second dimension) for each signal. unit_margin=1
because the shift amount is depend on the signal number.
An array with same dimensions as the input x
, but with
index shifted. The missing elements will be filled with NA
.
# Set ncores = 2 to comply to CRAN policy. Please don't run this line ravetools_threads(n_threads = 2L) x < matrix(1:10, nrow = 2, byrow = TRUE) z < shift_array(x, 2, 1, c(1,2)) y < NA * x y[1,1:4] = x[1,2:5] y[2,1:3] = x[2,3:5] # Check if z ang y are the same z  y # array case # x is Trial x Frequency x Time x < array(1:27, c(3,3,3)) # Shift time for each trial, amount is 1, 1, 0 shift_amount < c(1,1,0) z < shift_array(x, 3, 1, shift_amount) oldpar < par(mfrow = c(3, 2), mai = c(0.8, 0.6, 0.4, 0.1)) for( ii in 1:3 ){ image(t(x[ii, ,]), ylab = 'Frequency', xlab = 'Time', main = paste('Trial', ii)) image(t(z[ii, ,]), ylab = 'Frequency', xlab = 'Time', main = paste('Shifted amount:', shift_amount[ii])) } par(oldpar)
# Set ncores = 2 to comply to CRAN policy. Please don't run this line ravetools_threads(n_threads = 2L) x < matrix(1:10, nrow = 2, byrow = TRUE) z < shift_array(x, 2, 1, c(1,2)) y < NA * x y[1,1:4] = x[1,2:5] y[2,1:3] = x[2,3:5] # Check if z ang y are the same z  y # array case # x is Trial x Frequency x Time x < array(1:27, c(3,3,3)) # Shift time for each trial, amount is 1, 1, 0 shift_amount < c(1,1,0) z < shift_array(x, 3, 1, shift_amount) oldpar < par(mfrow = c(3, 2), mai = c(0.8, 0.6, 0.4, 0.1)) for( ii in 1:3 ){ image(t(x[ii, ,]), ylab = 'Frequency', xlab = 'Time', main = paste('Trial', ii)) image(t(z[ii, ,]), ylab = 'Frequency', xlab = 'Time', main = paste('Shifted amount:', shift_amount[ii])) } par(oldpar)
Create surface from 3Darray using marching cubes algorithm
vcg_isosurface( volume, threshold_lb = 0, threshold_ub = NA, vox_to_ras = diag(c(1, 1, 1, 1)) )
vcg_isosurface( volume, threshold_lb = 0, threshold_ub = NA, vox_to_ras = diag(c(1, 1, 1, 1)) )
volume 
a volume or a mask volume 
threshold_lb 
lowerbound threshold for creating the surface; default
is 
threshold_ub 
upperbound threshold for creating the surface; default
is 
vox_to_ras 
a 
A triangular mesh of class 'mesh3d'
if(is_not_cran()) { library(ravetools) data("left_hippocampus_mask") mesh < vcg_isosurface(left_hippocampus_mask) rgl_view({ rgl_call("mfrow3d", 1, 2) rgl_call("title3d", "Direct ISOSurface") rgl_call("shade3d", mesh, col = 2) rgl_call("next3d") rgl_call("title3d", "ISOSurface + Implicit Smooth") rgl_call("shade3d", vcg_smooth_implicit(mesh, degree = 2), col = 3) }) }
if(is_not_cran()) { library(ravetools) data("left_hippocampus_mask") mesh < vcg_isosurface(left_hippocampus_mask) rgl_view({ rgl_call("mfrow3d", 1, 2) rgl_call("title3d", "Direct ISOSurface") rgl_call("shade3d", mesh, col = 2) rgl_call("next3d") rgl_call("title3d", "ISOSurface + Implicit Smooth") rgl_call("shade3d", vcg_smooth_implicit(mesh, degree = 2), col = 3) }) }
Compute volume for manifold meshes
vcg_mesh_volume(mesh)
vcg_mesh_volume(mesh)
mesh 
triangular mesh of class 
The numeric volume of the mesh
# Initial mesh mesh < vcg_sphere() vcg_mesh_volume(mesh)
# Initial mesh mesh < vcg_sphere() vcg_mesh_volume(mesh)
Applies smoothing algorithms on a triangular mesh.
vcg_smooth_implicit( mesh, lambda = 0.2, use_mass_matrix = TRUE, fix_border = FALSE, use_cot_weight = FALSE, degree = 1L, laplacian_weight = 1 ) vcg_smooth_explicit( mesh, type = c("taubin", "laplace", "HClaplace", "fujiLaplace", "angWeight", "surfPreserveLaplace"), iteration = 10, lambda = 0.5, mu = 0.53, delta = 0.1 )
vcg_smooth_implicit( mesh, lambda = 0.2, use_mass_matrix = TRUE, fix_border = FALSE, use_cot_weight = FALSE, degree = 1L, laplacian_weight = 1 ) vcg_smooth_explicit( mesh, type = c("taubin", "laplace", "HClaplace", "fujiLaplace", "angWeight", "surfPreserveLaplace"), iteration = 10, lambda = 0.5, mu = 0.53, delta = 0.1 )
mesh 
triangular mesh stored as object of class 'mesh3d'. 
lambda 
In 
use_mass_matrix 
logical: whether to use mass matrix to keep the mesh
close to its original position (weighted per area distributed on vertices);
default is 
fix_border 
logical: whether to fix the border vertices of the mesh;
default is 
use_cot_weight 
logical: whether to use cotangent weight; default is

degree 
integer: degrees of 'Laplacian'; default is 
laplacian_weight 
numeric: weight when 
type 
method name of explicit smooth, choices are 
iteration 
number of iterations 
mu 
parameter for 
delta 
parameter for scaledependent 'Laplacian' smoothing or maximum allowed angle (in 'Radian') for deviation between surface preserving 'Laplacian'. 
An object of class "mesh3d" with:
vb 
vertex coordinates 
normals 
vertex normal vectors 
it 
triangular face index 
if(is_not_cran()) { # Prepare mesh with no normals data("left_hippocampus_mask") # Grow 2mm on each direction to fill holes volume < grow_volume(left_hippocampus_mask, 2) # Initial mesh mesh < vcg_isosurface(volume) # Start: examples rgl_view({ rgl_call("mfrow3d", 2, 4) rgl_call("title3d", "Naive ISOSurface") rgl_call("shade3d", mesh, col = 2) rgl_call("next3d") rgl_call("title3d", "Implicit Smooth") rgl_call("shade3d", col = 2, x = vcg_smooth_implicit(mesh, degree = 2)) rgl_call("next3d") rgl_call("title3d", "Explicit Smooth  taubin") rgl_call("shade3d", col = 2, x = vcg_smooth_explicit(mesh, "taubin")) rgl_call("next3d") rgl_call("title3d", "Explicit Smooth  laplace") rgl_call("shade3d", col = 2, x = vcg_smooth_explicit(mesh, "laplace")) rgl_call("next3d") rgl_call("title3d", "Explicit Smooth  angWeight") rgl_call("shade3d", col = 2, x = vcg_smooth_explicit(mesh, "angWeight")) rgl_call("next3d") rgl_call("title3d", "Explicit Smooth  HClaplace") rgl_call("shade3d", col = 2, x = vcg_smooth_explicit(mesh, "HClaplace")) rgl_call("next3d") rgl_call("title3d", "Explicit Smooth  fujiLaplace") rgl_call("shade3d", col = 2, x = vcg_smooth_explicit(mesh, "fujiLaplace")) rgl_call("next3d") rgl_call("title3d", "Explicit Smooth  surfPreserveLaplace") rgl_call("shade3d", col = 2, x = vcg_smooth_explicit(mesh, "surfPreserveLaplace")) }) }
if(is_not_cran()) { # Prepare mesh with no normals data("left_hippocampus_mask") # Grow 2mm on each direction to fill holes volume < grow_volume(left_hippocampus_mask, 2) # Initial mesh mesh < vcg_isosurface(volume) # Start: examples rgl_view({ rgl_call("mfrow3d", 2, 4) rgl_call("title3d", "Naive ISOSurface") rgl_call("shade3d", mesh, col = 2) rgl_call("next3d") rgl_call("title3d", "Implicit Smooth") rgl_call("shade3d", col = 2, x = vcg_smooth_implicit(mesh, degree = 2)) rgl_call("next3d") rgl_call("title3d", "Explicit Smooth  taubin") rgl_call("shade3d", col = 2, x = vcg_smooth_explicit(mesh, "taubin")) rgl_call("next3d") rgl_call("title3d", "Explicit Smooth  laplace") rgl_call("shade3d", col = 2, x = vcg_smooth_explicit(mesh, "laplace")) rgl_call("next3d") rgl_call("title3d", "Explicit Smooth  angWeight") rgl_call("shade3d", col = 2, x = vcg_smooth_explicit(mesh, "angWeight")) rgl_call("next3d") rgl_call("title3d", "Explicit Smooth  HClaplace") rgl_call("shade3d", col = 2, x = vcg_smooth_explicit(mesh, "HClaplace")) rgl_call("next3d") rgl_call("title3d", "Explicit Smooth  fujiLaplace") rgl_call("shade3d", col = 2, x = vcg_smooth_explicit(mesh, "fujiLaplace")) rgl_call("next3d") rgl_call("title3d", "Explicit Smooth  surfPreserveLaplace") rgl_call("shade3d", col = 2, x = vcg_smooth_explicit(mesh, "surfPreserveLaplace")) }) }
Simple 3dimensional sphere mesh
vcg_sphere(sub_division = 3L, normals = TRUE)
vcg_sphere(sub_division = 3L, normals = TRUE)
sub_division 
density of vertex in the resulting mesh 
normals 
whether the normal vectors should be calculated 
A 'mesh3d'
object
vcg_sphere()
vcg_sphere()
Sample a surface mesh uniformly
vcg_uniform_remesh( x, voxel_size = NULL, offset = 0, discretize = FALSE, multi_sample = FALSE, absolute_distance = FALSE, merge_clost = FALSE, verbose = TRUE )
vcg_uniform_remesh( x, voxel_size = NULL, offset = 0, discretize = FALSE, multi_sample = FALSE, absolute_distance = FALSE, merge_clost = FALSE, verbose = TRUE )
x 
surface 
voxel_size 
'voxel' size for space 'discretization' 
offset 
offset position shift of the new surface from the input 
discretize 
whether to use step function ( 
multi_sample 
whether to calculate multiple samples for more accurate
results (at the expense of more computing time) to remove artifacts; default
is 
absolute_distance 
whether an unsigned distance field should be
computed. When set to 
merge_clost 
whether to merge close vertices; default is 
verbose 
whether to verbose the progress; default is 
A triangular mesh of class 'mesh3d'
sphere < vcg_sphere() mesh < vcg_uniform_remesh(sphere, voxel_size = 0.45) if(is_not_cran()) { rgl_view({ rgl_call("mfrow3d", 1, 2) rgl_call("title3d", "Input") rgl_call("wire3d", sphere, col = 2) rgl_call("next3d") rgl_call("title3d", "Remeshed to 0.1mm edge distance") rgl_call("wire3d", mesh, col = 3) }) }
sphere < vcg_sphere() mesh < vcg_uniform_remesh(sphere, voxel_size = 0.45) if(is_not_cran()) { rgl_view({ rgl_call("mfrow3d", 1, 2) rgl_call("title3d", "Input") rgl_call("wire3d", sphere, col = 2) rgl_call("next3d") rgl_call("title3d", "Remeshed to 0.1mm edge distance") rgl_call("wire3d", mesh, col = 3) }) }
Update vertex normal
vcg_update_normals( mesh, weight = c("area", "angle"), pointcloud = c(10, 0), verbose = FALSE )
vcg_update_normals( mesh, weight = c("area", "angle"), pointcloud = c(10, 0), verbose = FALSE )
mesh 
triangular mesh or a pointcloud (matrix of 3 columns) 
weight 
method to compute pervertex normal vectors: 
pointcloud 
integer vector of length 2: containing optional parameters for normal calculation of point clouds; the first entry specifies the number of neighboring points to consider; the second entry specifies the amount of smoothing iterations to be performed. 
verbose 
whether to verbose the progress 
A 'mesh3d'
object with normal vectors.
if(is_not_cran()) { # Prepare mesh with no normal data("left_hippocampus_mask") mesh < vcg_isosurface(left_hippocampus_mask) mesh$normals < NULL # Start: examples new_mesh < vcg_update_normals(mesh, weight = "angle", pointcloud = c(10, 10)) rgl_view({ rgl_call("mfrow3d", 1, 2) rgl_call("shade3d", mesh, col = 2) rgl_call("next3d") rgl_call("shade3d", new_mesh, col = 2) }) }
if(is_not_cran()) { # Prepare mesh with no normal data("left_hippocampus_mask") mesh < vcg_isosurface(left_hippocampus_mask) mesh$normals < NULL # Start: examples new_mesh < vcg_update_normals(mesh, weight = "angle", pointcloud = c(10, 10)) rgl_view({ rgl_call("mfrow3d", 1, 2) rgl_call("shade3d", mesh, col = 2) rgl_call("next3d") rgl_call("shade3d", new_mesh, col = 2) }) }
Transform analog voltage signals with 'Morlet'
wavelets: complex wavelet kernels with $\pi/2$
phase
differences.
wavelet_kernels(freqs, srate, wave_num) morlet_wavelet( data, freqs, srate, wave_num, precision = c("float", "double"), trend = c("constant", "linear", "none"), signature = NULL, ... ) wavelet_cycles_suggest( freqs, frequency_range = c(2, 200), cycle_range = c(3, 20) )
wavelet_kernels(freqs, srate, wave_num) morlet_wavelet( data, freqs, srate, wave_num, precision = c("float", "double"), trend = c("constant", "linear", "none"), signature = NULL, ... ) wavelet_cycles_suggest( freqs, frequency_range = c(2, 200), cycle_range = c(3, 20) )
freqs 
frequency in which 
srate 
sample rate, number of time points per second 
wave_num 
desired number of cycles in wavelet kernels to balance the precision in time and amplitude (control the smoothness); positive integers are strongly suggested 
data 
numerical vector such as analog voltage signals 
precision 
the precision of computation; choices are

trend 
choices are 
signature 
signature to calculate kernel path to save, internally used 
... 
further passed to 
frequency_range 
frequency range to calculate, default is 2 to 200 
cycle_range 
number of cycles corresponding to 
wavelet_kernels
returns wavelet kernels to be
used for wavelet function; morlet_wavelet
returns a filebased array
if precision
is 'float'
, or a list of real and imaginary
arrays if precision
is 'double'
# generate sine waves time < seq(0, 3, by = 0.01) x < sin(time * 20*pi) + exp(time^2) * cos(time * 10*pi) plot(time, x, type = 'l') # freq from 1  15 Hz; wavelet using float precision freq < seq(1, 15, 0.2) coef < morlet_wavelet(x, freq, 100, c(2,3)) # to get coefficients in complex number from 110 time points coef[1:10, ] # power power < Mod(coef[])^2 # Power peaks at 5Hz and 10Hz at early stages # After 1.0 second, 5Hz component fade away image(power, x = time, y = freq, ylab = "frequency") # wavelet using double precision coef2 < morlet_wavelet(x, freq, 100, c(2,3), precision = "double") power2 < (coef2$real[])^2 + (coef2$imag[])^2 image(power2, x = time, y = freq, ylab = "frequency") # The maximum relative change of power with different precisions max(abs(power/power2  1)) # display kernels freq < seq(1, 15, 1) kern < wavelet_kernels(freq, 100, c(2,3)) print(kern) plot(kern)
# generate sine waves time < seq(0, 3, by = 0.01) x < sin(time * 20*pi) + exp(time^2) * cos(time * 10*pi) plot(time, x, type = 'l') # freq from 1  15 Hz; wavelet using float precision freq < seq(1, 15, 0.2) coef < morlet_wavelet(x, freq, 100, c(2,3)) # to get coefficients in complex number from 110 time points coef[1:10, ] # power power < Mod(coef[])^2 # Power peaks at 5Hz and 10Hz at early stages # After 1.0 second, 5Hz component fade away image(power, x = time, y = freq, ylab = "frequency") # wavelet using double precision coef2 < morlet_wavelet(x, freq, 100, c(2,3), precision = "double") power2 < (coef2$real[])^2 + (coef2$imag[])^2 image(power2, x = time, y = freq, ylab = "frequency") # The maximum relative change of power with different precisions max(abs(power/power2  1)) # display kernels freq < seq(1, 15, 1) kern < wavelet_kernels(freq, 100, c(2,3)) print(kern) plot(kern)