### UTILITIES FOR LOW-LEVEL MATH ###
#' Report time
#'
#' Provides a nicely formatted "estimated time left" in loops plus a summary
#' upon completion.
#' @param i current iteration
#' @param time_start time when the loop started running
#' @param nIter total number of iterations
#' @param reportEvery report progress every n iterations
#' @param jobs vector of length \code{nIter} specifying the relative difficulty
#' of each iteration. If not NULL, estimated time left takes into account
#' whether the jobs ahead will take more or less time than the jobs already
#' completed
#' @param prefix a string to print before "Done...", eg "Chain 1: "
#' @export
#' @examples
#' time_start = proc.time()
#' for (i in 1:100) {
#' Sys.sleep(i ^ 1.02 / 10000)
#' reportTime(i = i, time_start = time_start, nIter = 100,
#' jobs = (1:100) ^ 1.02, prefix = 'Chain 1: ')
#' }
#' \dontrun{
#' # Unknown number of iterations:
#' time_start = proc.time()
#' for (i in 1:20) {
#' Sys.sleep(i ^ 2 / 10000)
#' reportTime(i = i, time_start = time_start,
#' jobs = (1:20) ^ 2, reportEvery = 5)
#' }
#'
#' # when analyzing a bunch of audio files, their size is a good estimate
#' # of how long each will take to process
#' time_start = proc.time()
#' filenames = list.files('~/Downloads/temp', pattern = "*.wav|.mp3",
#' full.names = TRUE)
#' filesizes = file.info(filenames)$size
#' for (i in 1:length(filenames)) {
#' # ...do what you have to do with each file...
#' reportTime(i = i, nIter = length(filenames),
#' time_start = time_start, jobs = filesizes)
#' }
#' }
reportTime = function(
i,
time_start,
nIter = NULL,
reportEvery = NULL,
jobs = NULL,
prefix = ''
) {
time_diff = as.numeric((proc.time() - time_start)[3])
if (is.null(reportEvery))
reportEvery = ifelse(is.null(nIter),
1,
max(1, 10 ^ (floor(log10(nIter)) - 1)))
if (is.null(nIter)) {
# number of iter unknown, so we just report time elapsed
if (i %% reportEvery == 0) {
print(paste0(prefix, 'Completed ', i, ' iterations in ',
convert_sec_to_hms(time_diff)))
}
} else {
# we know how many iter, so we also report time left
if (i == nIter) {
time_total = convert_sec_to_hms(time_diff)
print(paste0(prefix, 'Completed ', i, ' iterations in ', time_total, '.'))
} else {
if (i %% reportEvery == 0 || i == 1) {
if (is.null(jobs)) {
# simply count iterations
time_left = time_diff / i * (nIter - i)
} else {
# take into account the expected time for each iteration
speed = time_diff / sum(jobs[1:i])
time_left = speed * sum(jobs[min((i + 1), nIter):nIter])
}
time_left_hms = convert_sec_to_hms(time_left)
print(paste0(prefix, 'Done ', i, ' / ', nIter, '; Estimated time left: ', time_left_hms))
}
}
}
}
#' Print time
#'
#' Internal soundgen function.
#'
#' Converts time in seconds to time in y m d h min s for pretty printing.
#' @param time_s time (s)
#' @param digits number of digits to preserve for s (1-60 s)
#' @return Returns a character string like "1 h 20 min 3 s"
#' @keywords internal
#' @examples
#' time_s = c(.0001, .01, .33, .8, 2.135, 5.4, 12, 250, 3721, 10000,
#' 150000, 365 * 24 * 3600 + 35 * 24 * 3600 + 3721)
#' soundgen:::convert_sec_to_hms(time_s)
#' soundgen:::convert_sec_to_hms(time_s, 2)
convert_sec_to_hms = function(time_s, digits = 0) {
if (!any(time_s > 1)) {
output = paste(round(time_s * 1000), 'ms')
} else {
len = length(time_s)
output = vector('character', len)
for (i in 1:len) {
# years = time_s %/% 31536000
# years_string = ifelse(years > 0, paste(years, 'y '), '')
#
# months = time_s %/% 2592000 - years * 12
# months_string = ifelse(months > 0, paste(months, 'm '), '')
days_string = hours_string = minutes_string = seconds_string = ms_string = ''
days = time_s[i] %/% 86400
if (days > 0) days_string = paste(days, 'd ')
hours = time_s[i] %/% 3600 - days * 24
if (hours > 0) hours_string = paste(hours, 'h ')
if (days == 0) {
minutes = time_s[i] %/% 60 - days * 1440 - hours * 60
if (minutes > 0) minutes_string = paste(minutes, 'min ')
if (hours == 0) {
seconds = time_s[i] - days * 86400 - hours * 3600 - minutes * 60
seconds_floor = floor(seconds)
if (seconds_floor > 0) seconds_string = paste(round(seconds, digits), 's ')
if (minutes == 0 & seconds_floor == 0) {
ms = (time_s[i] %% 1) * 1000
if (ms > 0) ms_string = paste(ms, 'ms')
}
}
}
output[i] = paste0(days_string, hours_string,
minutes_string, seconds_string, ms_string)
}
}
output = trimws(output)
return(output)
}
#' Convert Hz to semitones
#'
#' Converts from Hz to semitones above C-5 (~0.5109875 Hz) or another reference
#' frequency. This may not seem very useful, but note that this gives us a nice
#' logarithmic scale for generating natural pitch transitions.
#'
#' @seealso \code{\link{semitonesToHz}} \code{\link{HzToNotes}}
#'
#' @param h vector or matrix of frequencies (Hz)
#' @param ref frequency of the reference value (defaults to C-5, 0.51 Hz)
#' @export
#' @examples
#' s = HzToSemitones(c(440, 293, 115))
#' # to convert to musical notation
#' notesDict$note[1 + round(s)]
#' # note the "1 +": semitones ABOVE C-5, i.e. notesDict[1, ] is C-5
#'
#' # Any reference tone can be specified. For ex., for semitones above C0, use:
#' HzToSemitones(440, ref = 16.35)
#' # TIP: see notesDict for a table of Hz frequencies to musical notation
HzToSemitones = function(h, ref = 0.5109875) {
return(log2(h / ref) * 12)
}
#' Convert semitones to Hz
#'
#' Converts from semitones above C-5 (~0.5109875 Hz) or another reference
#' frequency to Hz. See \code{\link{HzToSemitones}}
#'
#' @seealso \code{\link{HzToSemitones}}
#'
#' @param s vector or matrix of frequencies (semitones above C0)
#' @param ref frequency of the reference value (defaults to C-5, 0.51 Hz)
#' @export
semitonesToHz = function(s, ref = 0.5109875) {
return(ref * 2 ^ (s / 12))
}
#' Convert Hz to notes
#'
#' Converts from Hz to musical notation like A4 - note A of the fourth octave
#' above C0 (16.35 Hz).
#'
#' @seealso \code{\link{notesToHz}} \code{\link{HzToSemitones}}
#'
#' @param h vector or matrix of frequencies (Hz)
#' @param showCents if TRUE, show cents to the nearest notes (cent = 1/100 of a
#' semitone)
#' @param A4 frequency of note A in the fourth octave (modern standard ISO 16 or
#' concert pitch = 440 Hz)
#' @export
#' @examples
#' HzToNotes(c(440, 293, 115, 16.35, 4))
#'
#' HzToNotes(c(440, 415, 80, 81), showCents = TRUE)
#' # 80 Hz is almost exactly midway (+49 cents) between D#2 and E2
#'
#' # Baroque tuning A415, half a semitone flat relative to concert pitch A440
#' HzToNotes(c(440, 415, 16.35), A4 = 415)
HzToNotes = function(h, showCents = FALSE, A4 = 440) {
ref = A4 / 861.0779
# C4 is 9 semitones below A4 (2^(9/12) * 512 = 861.0779), then another 9 octaves
# down to C-5 in notesDict
sem = log2(h / ref) * 12
round_sem = round(sem)
nearest_note = soundgen::notesDict$note[1 + round_sem]
if (showCents) {
cents = paste0(' ', as.character(round((sem - round_sem) * 100)))
cents[cents == ' 0'] = ''
nearest_note = paste0(nearest_note, cents)
}
return(nearest_note)
}
#' Convert notes to Hz
#'
#' Converts to Hz from musical notation like A4 - note A of the fourth octave
#' above C0 (16.35 Hz).
#'
#' @seealso \code{\link{HzToNotes}} \code{\link{HzToSemitones}}
#'
#' @param n vector or matrix of notes
#' @param A4 frequency of note A in the fourth octave (modern standard ISO 16 or
#' concert pitch = 440 Hz)
#' @export
#' @examples
#' notesToHz(c("A4", "D4", "A#2", "C0", "C-2"))
#'
#' # Baroque tuning A415, half a semitone flat relative to concert pitch A440
#' notesToHz(c("A4", "D4", "A#2", "C0", "C-2"), A4 = 415)
notesToHz = function(n, A4 = 440) {
soundgen::notesDict$freq[match(n, soundgen::notesDict$note)] * A4 / 440
}
#' Convert Hz to ERB rate
#'
#' Converts from Hz to the number of Equivalent Rectangular Bandwidths (ERBs)
#' below input frequency. See https://www2.ling.su.se/staff/hartmut/bark.htm and
#' https://en.wikipedia.org/wiki/Equivalent_rectangular_bandwidth
#'
#' @seealso \code{\link{ERBToHz}} \code{\link{HzToSemitones}}
#' \code{\link{HzToNotes}}
#'
#' @param h vector or matrix of frequencies (Hz)
#' @param method approximation to use
#' @export
#' @examples
#' HzToERB(c(-20, 20, 100, 440, 1000, NA))
#'
#' f = 20:20000
#' erb_lin = HzToERB(f, 'linear')
#' erb_quadratic = HzToERB(f, 'quadratic')
#' plot(f, erb_lin, log = 'x', type = 'l')
#' points(f, erb_quadratic, col = 'blue', type = 'l')
#'
#' # compare with the bark scale:
#' barks = tuneR::hz2bark(f)
#' points(f, barks / max(barks) * max(erb_lin),
#' col = 'red', type = 'l', lty = 2)
HzToERB = function(h,
method = c('linear', 'quadratic')[1]) {
if (method == 'linear') {
erb = 21.4 * log10(1 + .00437 * h)
} else if (method == 'quadratic') {
# erb = 11.17 * log( (h + 312) / (h + 14675)) + 43
erb = 11.17268 * log(1 + (46.06538 * h) / (h + 14678.49))
}
return(erb)
}
#' Convert Hz to ERB rate
#'
#' Converts from Hz to the number of Equivalent Rectangular Bandwidths (ERBs)
#' below input frequency. See https://www2.ling.su.se/staff/hartmut/bark.htm and
#' https://en.wikipedia.org/wiki/Equivalent_rectangular_bandwidth
#'
#' @seealso \code{\link{HzToERB}} \code{\link{HzToSemitones}}
#' \code{\link{HzToNotes}}
#'
#' @param e vector or matrix of frequencies in ERB rate
#' @param method approximation to use
#' @export
#' @examples
#' freqs_Hz = c(-20, 20, 100, 440, 1000, 20000, NA)
#' e_lin = HzToERB(freqs_Hz, 'linear')
#' ERBToHz(e_lin, 'linear')
#'
#' e_quad = HzToERB(freqs_Hz, 'quadratic')
#' ERBToHz(e_quad, 'quadratic')
ERBToHz = function(e,
method = c('linear', 'quadratic')[1]) {
if (method == 'linear') {
h = (10 ^ (e / 21.4) - 1) / .00437
} else if (method == 'quadratic') {
h = 676170.4 / (47.06538 - exp(e * 0.08950404)) - 14678.49
}
return(h)
}
#' Hz to mel
#'
#' Internal soundgen function: a temporary fix needed because tuneR::hz2mel
#' doesn't accept NAs.
#' @param f frequency, Hz
#' @param htk algrithm
#' @keywords internal
#' @examples
#' soundgen:::hz2mel(c(440, NA))
hz2mel = function(f, htk = FALSE) {
# if (!is.numeric(f) || f < 0)
# stop("frequencies have to be non-negative")
idx_nonFinite = which(!is.finite(f))
f[idx_nonFinite] = 0
if (htk) {
z <- 2595 * log10(1 + f/700)
}
else {
f_0 <- 0
f_sp <- 200/3
brkfrq <- 1000
brkpt <- (brkfrq - f_0)/f_sp
logstep <- exp(log(6.4)/27)
linpts <- (f < brkfrq)
z <- 0 * f
z[linpts] <- (f[linpts] - f_0)/f_sp
z[!linpts] <- brkpt + (log(f[!linpts]/brkfrq))/log(logstep)
}
z[idx_nonFinite] = NA
return(z)
}
#' Convert to dB
#'
#' Internal soundgen function.
#' @param x a vector of floats between 0 and 1 (exclusive, i.e. these are ratios)
#' @keywords internal
#' @examples
#' soundgen:::to_dB(c(.1, .5, .75, .9, .95, .99, .999, .9999))
to_dB = function(x) {
return(10 * log10(x / (1 - x)))
}
#' List depth
#'
#' Internal soundgen function
#'
#' Returns the depth of list structure. See https://stackoverflow.com/questions/13432863/determine-level-of-nesting-in-r
#' @param x any R object
#' @keywords internal
listDepth = function(x) ifelse(is.list(x), 1L + max(sapply(x, listDepth)), 0L)
#' Normalize 0 to 1
#'
#' Internal soundgen function
#'
#' Normalized input vector to range from 0 to 1
#' @param x numeric vector or matrix
#' @param na.rm if TRUE, removed NA's when calculating min/max for normalization
#' @param xmin,xmax min and max (to save time if already known)
#' @keywords internal
zeroOne = function(x, na.rm = FALSE, xmin = NULL, xmax = NULL) {
if (is.null(xmin)) xmin = min(x, na.rm = na.rm)
if (is.null(xmax)) xmax = max(x, na.rm = na.rm)
x = x - xmin
x = x / (xmax - xmin)
return(x)
}
#' log01
#'
#' Internal soundgen function
#'
#' Normalizes, log-transforms, and re-normalizes an input vector, so it ranges
#' from 0 to 1
#' @param v numeric vector
#' @keywords internal
#' @examples
#' v = exp(1:10)
#' soundgen:::log01(v)
log01 = function(v) {
v = v - min(v) + 1
v = log(v)
v = zeroOne(v)
return(v)
}
#' Entropy
#'
#' Returns Weiner or Shannon entropy of an input vector such as the spectrum of
#' a sound. Non-positive input values are converted to a small positive number
#' (convertNonPositive). If all elements are zero, returns NA.
#'
#' @param x vector of positive floats
#' @param type 'shannon' for Shannon (information) entropy, 'weiner' for Weiner
#' entropy
#' @param normalize if TRUE, Shannon entropy is normalized by the length of
#' input vector to range from 0 to 1. It has no affect on Weiner entropy
#' @param convertNonPositive all non-positive values are converted to
#' \code{convertNonPositive}
#' @export
#' @examples
#' # Here are four simplified power spectra, each with 9 frequency bins:
#' s = list(
#' c(rep(0, 4), 1, rep(0, 4)), # a single peak in spectrum
#' c(0, 0, 1, 0, 0, .75, 0, 0, .5), # perfectly periodic, with 3 harmonics
#' rep(0, 9), # a silent frame
#' rep(1, 9) # white noise
#' )
#'
#' # Weiner entropy is ~0 for periodic, NA for silent, 1 for white noise
#' sapply(s, function(x) round(getEntropy(x), 2))
#'
#' # Shannon entropy is ~0 for periodic with a single harmonic, moderate for
#' # periodic with multiple harmonics, NA for silent, highest for white noise
#' sapply(s, function(x) round(getEntropy(x, type = 'shannon'), 2))
#'
#' # Normalized Shannon entropy - same but forced to be 0 to 1
#' sapply(s, function(x) round(getEntropy(x,
#' type = 'shannon', normalize = TRUE), 2))
getEntropy = function(x,
type = c('weiner', 'shannon')[1],
normalize = FALSE,
convertNonPositive = 1e-10) {
if (sum(x) == 0) return (NA) # empty frames
x = ifelse(x <= 0, convertNonPositive, x) # otherwise log0 gives NaN
if (type == 'weiner') {
geom_mean = exp(mean(log(x)))
ar_mean = mean(x)
entropy = geom_mean / ar_mean
} else if (type == 'shannon') {
p = x / sum(x)
if (normalize) {
entropy = -sum(p * log2(p) / log(length(p)) * log(2))
} else { # unnormalized
entropy = -sum(p * log2(p))
}
} else {
stop('Implemented entropy types: "shannon" or "weiner"')
}
return (entropy)
}
#' Random draw from a truncated normal distribution
#'
#' \code{rnorm_truncated} generates random numbers from a normal distribution
#' using rnorm(), but forced to remain within the specified low/high bounds. All
#' proposals outside the boundaries (exclusive) are discarded, and the sampling
#' is repeated until there are enough values within the specified range. Fully
#' vectorized.
#'
#' @param n the number of values to return
#' @param mean the mean of the normal distribution from which values are
#' generated (vector of length 1 or n)
#' @param sd the standard deviation of the normal distribution from which values
#' are generated (vector of length 1 or n)
#' @param low,high exclusive lower and upper bounds ((vectors of length 1 or n))
#' @param roundToInteger boolean vector of length 1 or n. If TRUE, the
#' corresponding value is rounded to the nearest integer.
#' @inheritParams soundgen
#' @return A vector of length n.
#' @keywords internal
#' @examples
#' soundgen:::rnorm_truncated(n = 3, mean = 10, sd = 5, low = 7, high = NULL,
#' roundToInteger = c(TRUE, FALSE, FALSE))
#' soundgen:::rnorm_truncated(n = 9, mean = c(10, 50, 100), sd = c(5, 0, 20),
#' roundToInteger = TRUE) # vectorized
#' # in case of conflicts between mean and bounds, either adjust the mean:
#' soundgen:::rnorm_truncated(n = 3, mean = 10, sd = .1,
#' low = c(15, 0, 0), high = c(100, 100, 8), invalidArgAction = 'adjust')
#' #... or ignore the boundaries
#' soundgen:::rnorm_truncated(n = 3, mean = 10, sd = .1,
#' low = c(15, 0, 0), high = c(100, 100, 8), invalidArgAction = 'ignore')
rnorm_truncated = function(n = 1,
mean = 0,
sd = 1,
low = NULL,
high = NULL,
roundToInteger = FALSE,
invalidArgAction = c('adjust', 'abort', 'ignore')[1]) {
if (length(mean) < n) mean = spline(mean, n = n)$y
if (length(sd) < n) sd = spline(sd, n = n)$y
if (any(!is.finite(sd))) sd = mean / 10
sd[sd < 0] = 0
if (sum(sd != 0) == 0) {
out = mean
out[roundToInteger] = round(out[roundToInteger], 0)
return (out)
}
if (is.null(low) & is.null(high)) {
out = rnorm(n, mean, sd)
out[roundToInteger] = round (out[roundToInteger], 0)
return (out)
}
if (is.null(low)) low = rep(-Inf, n)
if (is.null(high)) high = rep(Inf, n)
if (length(low) == 1) low = rep(low, n)
if (length(high) == 1) high = rep(high, n)
if (any(mean > high | mean < low)) {
warning(paste('Some of the specified means are outside the low/high bounds!',
'Mean =', paste(head(mean, 3), collapse = ', '),
'Low =', paste(head(low, 3), collapse = ', '),
'High = ', paste(head(high, 3), collapse = ', ')))
if (invalidArgAction == 'abort') {
stop('Aborting rnorm_truncated()')
} else if (invalidArgAction == 'adjust') {
mean[mean < low] = low[mean < low]
mean[mean > high] = high[mean > high]
} else if (invalidArgAction == 'ignore') {
low = rep(-Inf, n)
high = rep(Inf, n)
}
}
out = rnorm(n, mean, sd)
out[roundToInteger] = round(out[roundToInteger], 0)
for (i in 1:n) {
while (out[i] < low[i] | out[i] > high[i]) {
out[i] = rnorm(1, mean[i], sd[i]) # repeat until a suitable value is generated
out[roundToInteger] = round(out[roundToInteger], 0)
}
}
return(out)
}
#' Modified mode
#'
#' Internal soundgen function
#'
#' Internal helper function for spectral (~BaNa) pitch tracker. NOT quite the same as simply mode(x).
#' @param x numeric vector
#' @keywords internal
#' @examples
#' soundgen:::Mode(c(1, 2, 3)) # if every element is unique, return the smallest
#' soundgen:::Mode(c(1, 2, 2, 3))
Mode = function(x) {
x = sort(x)
ux = unique(x)
if (length(ux) < length(x)) {
return(ux[which.max(tabulate(match(x, ux)))])
} else {
# if every element is unique, return the smallest
return(x[1])
}
}
#' Random walk
#'
#' Generates a random walk with flexible control over its range, trend, and
#' smoothness. It works by calling stats::rnorm at each step and taking a
#' cumulative sum of the generated values. Smoothness is controlled by initially
#' generating a shorter random walk and upsampling.
#' @param len an integer specifying the required length of random walk. If len
#' is 1, returns a single draw from a gamma distribution with mean=1 and
#' sd=rw_range
#' @param rw_range the upper bound of the generated random walk (the lower bound
#' is set to 0)
#' @param rw_smoothing specifies the amount of smoothing, basically the number
#' of points used to construct the rw as a proportion of len, from 0 (no
#' smoothing) to 1 (maximum smoothing to a straight line)
#' @param method specifies the method of smoothing: either linear interpolation
#' ('linear', see stats::approx) or cubic splines ('spline', see
#' stats::spline)
#' @param trend mean of generated normal distribution (vectors are also
#' acceptable, as long as their length is an integer multiple of len). If
#' positive, the random walk has an overall upwards trend (good values are
#' between 0 and 0.5 or -0.5). Trend = c(1,-1) gives a roughly bell-shaped rw
#' with an upward and a downward curve. Larger absolute values of trend
#' produce less and less random behavior
#' @return Returns a numeric vector of length len and range from 0 to rw_range.
#' @export
#' @examples
#' plot(getRandomWalk(len = 1000, rw_range = 5, rw_smoothing = 0))
#' plot(getRandomWalk(len = 1000, rw_range = 5, rw_smoothing = .2))
#' plot(getRandomWalk(len = 1000, rw_range = 5, rw_smoothing = .95))
#' plot(getRandomWalk(len = 1000, rw_range = 5, rw_smoothing = .99))
#' plot(getRandomWalk(len = 1000, rw_range = 5, rw_smoothing = 1))
#' plot(getRandomWalk(len = 1000, rw_range = 15,
#' rw_smoothing = .2, trend = c(.1, -.1)))
#' plot(getRandomWalk(len = 1000, rw_range = 15,
#' rw_smoothing = .2, trend = c(15, -1)))
getRandomWalk = function(len,
rw_range = 1,
rw_smoothing = .2,
method = c('linear', 'spline')[2],
trend = 0) {
if (len < 2)
return(rgamma(1, 1 / rw_range ^ 2, 1 / rw_range ^ 2))
# generate a random walk (rw) of length depending on rw_smoothing,
# then linear extrapolation to len
# n = floor(max(2, 2 ^ (1 / rw_smoothing)))
# n = 2 + (len - 1) / (1 + exp(10 * (rw_smoothing - .1)))
n = max(2, len * (1 - rw_smoothing))
if (length(trend) > 1) {
n = round(n / 2, 0) * 2 # force to be even
trend_short = rep(trend, each = n / length(trend))
# for this to work, length(trend) must be a multiple of n.
# In practice, specify trend of length 2
} else {
trend_short = trend
}
if (n > len) {
rw_long = cumsum(rnorm(len, trend_short)) # just a rw of length /len/
} else {
# get a shorter sequence and extrapolate, thus achieving
# more or less smoothing
rw_short = cumsum(rnorm(n, trend_short)) # plot(rw_short, type = 'l')
if (method == 'linear') {
rw_long = approx(rw_short, n = len)$y
} else if (method == 'spline') {
rw_long = spline(rw_short, n = len)$y
}
} # plot (rw_long, type = 'l')
# normalize
rw_normalized = rw_long - min(rw_long)
rw_normalized = rw_normalized / max(abs(rw_normalized)) * rw_range
return(rw_normalized)
}
#' Discrete random walk
#'
#' Takes a continuous random walk and converts it to continuous epochs of
#' repeated values 0/1/2, each at least minLength points long. 0/1/2 correspond
#' to different noise regimes: 0 = no noise, 1 = subharmonics, 2 = subharmonics
#' and jitter/shimmer.
#' @param rw a random walk generated by \code{\link{getRandomWalk}} (expected
#' range 0 to 100)
#' @param nonlinBalance a number between 0 to 100: 0 = returns all zeros;
#' 100 = returns all twos
#' @param minLength the mimimum length of each epoch
#' @param q1,q2 cutoff points for transitioning from regime 0 to 1 (q1) or from
#' regime 1 to 2 (q2). See noiseThresholdsDict for defaults
#' @param plot if TRUE, plots the random walk underlying nonlinear regimes
#' @return Returns a vector of integers (0/1/2) of the same length as rw.
#' @export
#' @examples
#' rw = getRandomWalk(len = 100, rw_range = 100, rw_smoothing = .2)
#' r = getIntegerRandomWalk(rw, nonlinBalance = 75,
#' minLength = 10, plot = TRUE)
#' r = getIntegerRandomWalk(rw, nonlinBalance = 15,
#' q1 = 30, q2 = 70,
#' minLength = 10, plot = TRUE)
getIntegerRandomWalk = function(rw,
nonlinBalance = 50,
minLength = 50,
q1 = NULL,
q2 = NULL,
plot = FALSE) {
len = length(rw)
# calculate thresholds for different noise regimes
if (is.null(q1)) {
q1 = noiseThresholdsDict$q1[nonlinBalance + 1]
# +1 b/c the rows indices in noiseThresholdsDict start from 0, not 1
}
if (is.null(q2)) {
q2 = noiseThresholdsDict$q2[nonlinBalance + 1]
}
if (nonlinBalance == 0) {
rw_bin = rep(0, len)
} else if (nonlinBalance == 100) {
rw_bin = rep(2, len)
} else {
# convert continuous rw to discrete epochs based on q1 and q2 thresholds
rw_bin = rep(0, len)
rw_bin[which(rw > q1)] = 1
rw_bin[which(rw > q2)] = 2 # plot (rw_bin, ylim=c(0,2))
# make sure each epoch is long enough
rw_bin = clumper(rw_bin, minLength = minLength)
}
if (plot) {
rw_bin_100 = rw_bin
rw_bin_100[rw_bin_100 == 1] = q1
rw_bin_100[rw_bin_100 == 2] = q2
plot(rw, ylim = c(0, 110), type = 'l', lwd = 1,
xlab = 'Time', ylab = 'Latent non-linearity', main = 'Random walk')
points(rw_bin_100, type = 'l', lwd = 4, col = 'blue')
lines(x = c(0, 100), y = c(q1, q1), lty = 3, lwd = 2, col = 'red')
text(x = 0, y = q1 + 2, labels = 'subh', pos = 4)
lines(x = c(0, 100), y = c(q2, q2), lty = 3, lwd = 2, col = 'red')
text(x = 0, y = q2 + 2, labels = 'subh + jitter', pos = 4)
}
return(rw_bin)
}
#' Resize vector to required length
#'
#' Internal soundgen function.
#'
#' Adjusts a vector to match the required length by either trimming one or both
#' ends or padding them with zeros.
#' @param myseq input vector
#' @param len target length
#' @param padDir specifies the affected side. For padding, it is the side on
#' which new elements will be added. For trimming, this is the side that will
#' be trimmed. Defaults to 'central'
#' @param padWith if the vector needs to be padded to match the required length,
#' what should it be padded with? Defaults to 0
#' @return Returns the modified vector of the required length.
#' @keywords internal
#' @examples
#' soundgen:::matchLengths(c(1, 2, 3), len = 5)
#' soundgen:::matchLengths(3:7, len = 3)
#' # trimmed on the left
#' soundgen:::matchLengths(3:7, len = 3, padDir = 'left')
#' # padded with zeros on the left
#' soundgen:::matchLengths(3:7, len = 10, padDir = 'left')
#' #' # trimmed on the right
#' soundgen:::matchLengths(3:7, len = 3, padDir = 'right')
#' # padded with zeros on the right
#' soundgen:::matchLengths(3:7, len = 10, padDir = 'right')
matchLengths = function(myseq,
len,
padDir = c('left', 'right', 'central')[3],
padWith = 0) {
# padDir specifies where to cut/add zeros ('left' / 'right' / 'central')
if (length(myseq) == len) return(myseq)
if (padDir == 'central') {
if (length(myseq) < len) {
myseq = c(rep(padWith, len), myseq, rep(padWith, len))
# for padding, first add a whole lot of zeros and then trim using the same
# algorithm as for trimming
}
halflen = len / 2
center = (1 + length(myseq)) / 2
start = ceiling(center - halflen)
myseq = myseq[start:(start + len - 1)]
} else if (padDir == 'left') {
if (length(myseq) > len) {
myseq = myseq [(length(myseq) - len + 1):length(myseq)]
} else {
myseq = c(rep(padWith, (len - length(myseq))), myseq)
}
} else if (padDir == 'right') {
if (length(myseq) > len) {
myseq = myseq [1:len]
} else {
myseq = c(myseq, rep(padWith, (len - length(myseq))))
}
}
return (myseq)
}
#' Match number of columns
#'
#' Internal soundgen function
#'
#' Adds columns of new values (eg zeros or NAs) to a matrix, so that the new
#' number of columns = \code{len}
#' @inheritParams matchLengths
#' @param matrix_short input matrix
#' @param nCol the required number of columns
#' @param padWith the value to pad with, normally \code{0} or \code{NA}
#' @keywords internal
#' @examples
#' a = matrix(1:9, nrow = 3)
#' soundgen:::matchColumns(a, nCol = 6, padWith = NA, padDir = 'central')
#' soundgen:::matchColumns(a, nCol = 6, padWith = 0, padDir = 'central')
#' soundgen:::matchColumns(a, nCol = 6, padWith = NA, padDir = 'left')
#' soundgen:::matchColumns(a, nCol = 6, padWith = 'a', padDir = 'right')
#' soundgen:::matchColumns(a, nCol = 2)
matchColumns = function (matrix_short,
nCol,
padWith = 0,
padDir = 'central',
interpol = c("approx", "spline")[1]) {
if (ncol(matrix_short) > nCol) {
# downsample
new = interpolMatrix(matrix_short, nc = nCol, interpol = interpol)
} else {
# pad with 0 / NA / etc
if (is.null(colnames(matrix_short))) {
col_short = 1:ncol(matrix_short)
} else {
col_short = colnames(matrix_short)
}
col_long = matchLengths(col_short, nCol,
padDir = padDir, padWith = NA)
new = matrix(padWith,
nrow = nrow(matrix_short),
ncol = length(col_long))
colnames(new) = col_long
# paste the old matrix where it belongs
new[, !is.na(colnames(new))] = matrix_short
}
return(new)
}
#' Add overlapping vectors
#'
#' Adds two partly overlapping vectors, such as two waveforms, to produce a
#' longer vector. The location at which vector 2 is pasted is defined by
#' insertionPoint. Algorithm: both vectors are padded with zeros to match in
#' length and then added. All NA's are converted to 0.
#'
#' @seealso \code{\link{soundgen}}
#'
#' @param v1,v2 numeric vectors
#' @param insertionPoint the index of element in vector 1 at which vector 2 will
#' be inserted (any integer, can also be negative)
#' @param normalize if TRUE, the output is normalized to range from -1 to +1
#' @export
#' @examples
#' v1 = 1:6
#' v2 = rep(100, 3)
#' addVectors(v1, v2, insertionPoint = 5, normalize = FALSE)
#' addVectors(v1, v2, insertionPoint = -4, normalize = FALSE)
#' # note the asymmetry: insertionPoint refers to the first arg
#' addVectors(v2, v1, insertionPoint = -4, normalize = FALSE)
#'
#' v3 = rep(100, 15)
#' addVectors(v1, v3, insertionPoint = -4, normalize = FALSE)
#' addVectors(v2, v3, insertionPoint = 7, normalize = FALSE)
addVectors = function(v1, v2, insertionPoint = 1, normalize = TRUE) {
if (!is.numeric(v1)) stop(paste('Non-numeric v1:', head(v1)))
if (!is.numeric(v2)) stop(paste('Non-numeric v2:', head(v2)))
v1[is.na(v1)] = 0
v2[is.na(v2)] = 0
# align left ends
if (insertionPoint > 1) {
pad_left = insertionPoint - 1
v2 = c(rep(0, insertionPoint), v2)
} else if (insertionPoint < 1) {
pad_left = 1 - insertionPoint
v1 = c(rep(0, pad_left), v1)
}
# equalize lengths
l1 = length(v1)
l2 = length(v2)
len_dif = l2 - l1
if (len_dif > 0) {
v1 = c(v1, rep(0, len_dif))
} else if (len_dif < 0) {
v2 = c(v2, rep(0, -len_dif))
}
# add and normalize
out = v1 + v2
if (normalize) out = out / max(abs(out))
return(out)
}
#' Clump a sequence into large segments
#'
#' Internal soundgen function.
#'
#' \code{clumper} makes sure each homogeneous segment in a sequence is at least
#' minLength long. Called by getIntegerRandomWalk(), getVocalFry(),
#' naiveBayes(), etc. Algorithm: find the epochs shorter than minLength, merge
#' max 1/4 of them with the largest neighbor, and repeat recursively until all
#' epochs are at least minLength long. minLength can be a vector, in which case
#' it is assumed to change over time.
#' @keywords internal
#' @param x a vector: anything that can be converted into an integer to call
#' diff(): factors, integers, characters, booleans
#' @param minLength the minimum length of a segment (interger or vector)
#' @return Returns the original sequence x transformed to homogeneous segments
#' of required length, with the original class (e.g. character or factor).
#' @keywords internal
#' @examples
#' s = c(1,3,2,2,2,0,0,4,4,1,1,1,1,1,3,3)
#' soundgen:::clumper(s, 2)
#' soundgen:::clumper(s, 3)
#' soundgen:::clumper(1:5, 10)
#' soundgen:::clumper(c('a','a','a','b','b','c','c','c','a','c'), 3)
#' soundgen:::clumper(x = c(1,2,1,2,1,1,1,1,3,1), minLength = c(1, 1, 1, 3))
#' soundgen:::clumper(as.factor(c('A','B','B','C')), 2)
clumper = function(x, minLength, n = length(x)) {
mode_orig = mode(x)
minLength = round(minLength) # just in case it's not all integers
if (max(minLength) < 2) return(x)
n = length(x)
if (n <= minLength[1]) {
# just repeat the most common element
tbx = table(x)
mode_x = names(tbx)[which.max(tbx)]
x = rep(mode_x, n)
mode(x) = mode_orig
return(x)
}
# find indices of changes in the composition of x (epoch boundaries)
if (inherits(x, 'character')) {
x_int = as.integer(as.factor(x))
} else {
x_int = as.integer(x)
}
idx_change = which(diff(as.numeric(x_int)) != 0)
if (length(idx_change) < 1) {
return(x)
}
dur_epochs = diff(unique(c(0, idx_change, n)))
n_epochs = length(dur_epochs)
if (length(minLength) == 1) {
minLength_ups = minLength
} else {
minLength_ups = round(approx(minLength, n = n_epochs)$y)
}
short_epochs = which(dur_epochs < minLength_ups)
short_epochs = short_epochs[order(dur_epochs[short_epochs])]
nse = length(short_epochs)
# limit the number of epochs that are to be modified per iteration - here set
# to max 1/4 of the total n_epochs (starting with the shortest - this is
# crucial!)
n_max = max(2, ceiling(n_epochs / 4))
if (nse > n_max) {
short_epochs = order(dur_epochs)[1:n_max]
nse = n_max
}
if (nse > 0) {
idx = c(1, idx_change + 1)
for (e in short_epochs) {
# calculate the indices within vector x corresponding to this short epoch
if (e < n_epochs) {
idx_epoch = idx[e]:(idx[e + 1] - 1)
} else {
idx_epoch = idx[e]:n
}
# merge the short epoch with the longer neighbor (ie replace original
# elements with those of the longest adjacent epoch)
dur_left = ifelse(e > 1, dur_epochs[e - 1], 0)
dur_right = ifelse(e < n_epochs, dur_epochs[e + 1], 0)
if (dur_left >= dur_right) {
x[idx_epoch] = x[idx_epoch[1] - 1]
} else if (dur_left < dur_right) {
x[idx_epoch] = x[tail(idx_epoch, 1) + 1]
}
}
# call clumper() recursively until no short epochs are left
x = clumper(x, minLength, n = n) # just to avoid recalculating n every time
}
return(x)
}
#' Simple peak detection
#'
#' Internal soundgen function.
#'
#' @param x input vector
#' @param threshold threshold for peak detection
#' @keywords internal
#' @examples
#' soundgen:::isCentral.localMax(c(1,1,3,2,1), 2.5)
isCentral.localMax = function(x, threshold) {
middle = ceiling(length(x) / 2)
return(x[middle] > threshold && which.max(x) == middle)
}
#' Get sigmoid filter
#'
#' Internal soundgen function.
#'
#' Produces a filter for amplitude modulation ranging from clicks to
#' approximately a sine wave to reversed clicks (small episodes of silence). The
#' filter is made from concatenated sigmoids and their mirror reflections.
#' @return Returns a vector of length \code{len} and range from 0 to 1
#' @param len the length of output vector
#' @param samplingRate the sampling rate of the output vector, Hz
#' @param freq the frequency of amplitude modulation, Hz (numeric vector)
#' @param shape 0 = ~sine, -1 = clicks, +1 = notches (NB: vice versa in
#' soundgen!); numeric vector of length 1 or the same length as \code{freq}
#' @param spikiness amplifies the effect of the "shape" parameter;
#' numeric vector of length 1 or the same length as \code{freq}
#' @keywords internal
#' @examples
#' for (shape in c(0, -.1, .1, -1, 1)) {
#' s = soundgen:::getSigmoid(shape = shape, len = 1000, samplingRate = 500, freq = 2)
#' plot(s, type = 'l', main = paste('shape =', shape), xlab = '', ylab = '')
#' }
getSigmoid = function(len,
samplingRate = 16000,
freq = 5,
shape = 0,
spikiness = 1) {
# print(c(len, freq))
if (length(freq) > 1 | length(shape) > 1 | length(spikiness) > 1) {
# get preliminary frequency contour to estimate how many cycles are needed
if (length(freq) != len) {
freqContour_prelim = getSmoothContour(
anchors = freq,
len = 100,
valueFloor = 0.001,
method = 'spline'
)
} else {
freqContour_prelim = freq
}
# plot(freqContour_prelim, type = 'l')
n = ceiling(len / samplingRate / mean(1 / freqContour_prelim))
# get actual contours
freqContour = getSmoothContour(anchors = freq, len = n,
valueFloor = 0.001, method = 'spline')
shapeContour = getSmoothContour(anchors = shape, len = n,
method = 'spline')
spikinessContour = getSmoothContour(anchors = spikiness,
len = n, method = 'spline')
# set up par vectors
from = -exp(-shapeContour * spikinessContour)
to = exp(shapeContour * spikinessContour)
slope = exp(abs(shapeContour)) * 5 # close to sine
# create one cycle at a time
out = vector()
i = 1
while (length(out) < len) {
a = seq(from = from[i], to = to[i],
length.out = samplingRate / freqContour[i] / 2)
b = 1 / (1 + exp(-a * slope[i]))
b = zeroOne(b) # plot(b, type = 'l')
out = c(out, b, rev(b))
i = ifelse(i + 1 > n, n, i + 1) # repeat the last value if we run out of cycles
}
out = out[1:len]
} else {
# if pars are static, we can take a shortcut (~50 times faster)
from = -exp(-shape * spikiness)
to = exp(shape * spikiness)
slope = exp(abs(shape)) * 5 # close to sine
a = seq(from = from, to = to, length.out = samplingRate / freq / 2)
b = 1 / (1 + exp(-a * slope))
b = zeroOne(b) # plot(b, type = 'l')
out = rep(c(b, rev(b)), length.out = len)
}
# plot(out, type = 'l')
return(out)
}
#' Report CI
#'
#' Internal soundgen function
#'
#' Takes a numeric vector or matrix with three elements / columns: fit, lower
#' quantile from a CI, and upper quantile from a CI. For each row, it prints the
#' result as fit and CI in square brackets
#' @param n numeric vector or matrix
#' @param digits number of decimal points to preserve
#' @keywords internal
#' @examples
#' n = rnorm(100)
#' soundgen:::reportCI(quantile(n, probs = c(.5, .025, .975)))
#'
#' a = data.frame(fit = c(3, 5, 7),
#' lwr = c(1, 4, 6.5),
#' upr = c(5, 6, 7.1))
#' soundgen:::reportCI(a, 1)
#' soundgen:::reportCI(a, 1, ' cm')
#' soundgen:::reportCI(a, 1, '%, 95% CI')
reportCI = function(n, digits = 2, suffix = NULL) {
if (is.data.frame(n)) n = as.matrix(n)
n = round(n, digits)
if (is.matrix(n)) {
out = matrix(NA, nrow = nrow(n))
rownames(out) = rownames(n)
for (i in 1:nrow(n)) {
out[i, ] = reportCI(n[i, ], digits = digits, suffix = suffix)
}
out
} else {
paste0(n[1], suffix, ' [', n[2], ', ', n[3], ']')
}
}
#' Interpolate matrix
#'
#' Internal soundgen function
#'
#' Performs a chosen type of interpolation (linear or smooth) across both rows
#' and columns of a matrix, in effect up- or downsampling a matrix to required
#' dimensions. Rownames and colnames are also interpolated as needed.
#' @param m input matrix of numeric values
#' @param nr,nc target dimensions
#' @param interpol interpolation method ('approx' for linear, 'spline' for
#' spline)
#' @keywords internal
#' @examples
#' m = matrix(1:12 + rnorm(12, 0, .2), nrow = 3)
#' rownames(m) = 1:3; colnames(m) = 1:4
#' soundgen:::interpolMatrix(m) # just returns the original
#' soundgen:::interpolMatrix(m, nr = 10, nc = 7)
#' soundgen:::interpolMatrix(m, nr = 10, nc = 7, interpol = 'spline')
#' soundgen:::interpolMatrix(m, nr = 2, nc = 7)
#' soundgen:::interpolMatrix(m, nr = 2, nc = 3)
#'
#' # input matrices can have a single row/column
#' soundgen:::interpolMatrix(matrix(1:5, nrow = 1), nc = 9)
#' soundgen:::interpolMatrix(matrix(1:5, ncol = 1), nr = 5, nc = 3)
interpolMatrix = function(m,
nr = NULL,
nc = NULL,
interpol = c("approx", "spline")[1]) {
if (!is.matrix(m)) {
m = as.matrix(m)
warning('non-matrix input m: converting to matrix')
}
nr0 = nrow(m)
nc0 = ncol(m)
if (is.null(nr)) nr = nr0
if (is.null(nc)) nc = nc0
if (nr == nr0 & nc == nc0) return(m)
# if (nr < 2) stop('nr must be >1')
# if (nc < 2) stop('nc must be >1')
isComplex = is.complex(m[1, 1])
# # Downsample rows if necessary
# if (nr0 > nr) {
# m = m[seq(1, nr0, length.out = nr),, drop = FALSE]
# }
#
# # Downsample columns if necessary
# if (nc0 > nc) {
# m = m[, seq(1, nc0, length.out = nc), drop = FALSE]
# }
# Interpolate rows if necessary
if (nr0 != nr) {
if (nr0 == 1) {
temp = matrix(rep(m, nr), nrow = nr, byrow = TRUE)
} else {
temp = matrix(1, nrow = nr, ncol = nc0)
for (c in 1:nc0) {
if (isComplex) {
# approx doesn't work with complex numbers properly, so we treat the
# Re and Im parts separately
temp_re = approx(x = Re(m[, c]), n = nr)$y
temp_im = approx(x = Im(m[, c]), n = nr)$y
temp[, c] = complex(real = temp_re, imaginary = temp_im)
} else {
temp[, c] = do.call(interpol, list(x = m[, c], n = nr))$y
}
}
}
if (!is.null(rownames(m))) {
rnms = as.numeric(rownames(m))
rownames(temp) = do.call(interpol, list(x = rnms, n = nr))$y
}
} else {
temp = m
rownames(temp) = rownames(m)
}
colnames(temp) = colnames(m)
# Interpolate columns if necessary
if (nc0 != nc) {
if (nc0 == 1) {
out = matrix(rep(temp[, 1], nc), ncol = nc, byrow = FALSE)
} else {
out = matrix(1, nrow = nr, ncol = nc)
for (r in 1:nr) {
if (isComplex) {
temp_re = approx(x = Re(temp[r, ]), n = nc)$y
temp_im = approx(x = Im(temp[r, ]), n = nc)$y
out[r, ] = complex(real = temp_re, imaginary = temp_im)
} else {
out[r, ] = do.call(interpol, list(x = temp[r, ], n = nc))$y
}
}
}
if (!is.null(colnames(m))) {
cnms = as.numeric(colnames(m))
try_colnames = try(do.call(interpol, list(x = cnms, n = nc))$y,
silent = TRUE)
if (!inherits(try_colnames, 'try-error')) {
colnames(out) = try_colnames
}
}
} else {
out = temp
colnames(out) = colnames(temp)
}
rownames(out) = rownames(temp)
return(out)
}
#' sampleModif
#'
#' Internal soundgen function
#'
#' Same as \code{\link[base]{sample}}, but without defaulting to x = 1:x if
#' length(x) = 1. See
#' https://stackoverflow.com/questions/7547758/using-sample-with-sample-space-size-1
#' @param x vector
#' @param ... other arguments passed to \code{sample}
#' @keywords internal
#' @examples
#' soundgen:::sampleModif(x = 3, n = 1)
#' # never returns 1 or 2: cf. sample(x = 3, n = 1)
sampleModif = function(x, ...) x[sample.int(length(x), ...)]
#' Gaussian smoothing in 2D
#'
#' Takes a matrix of numeric values and smoothes it by convolution with a
#' symmetric Gaussian window function.
#'
#' @seealso \code{\link{modulationSpectrum}}
#'
#' @return Returns a numeric matrix of the same dimensions as input.
#' @param m input matrix (numeric, on any scale, doesn't have to be square)
#' @param kernelSize the size of the Gaussian kernel, in points
#' @param kernelSD the SD of the Gaussian kernel relative to its size (.5 = the
#' edge is two SD's away)
#' @param action 'blur' = kernel-weighted average, 'unblur' = subtract
#' kernel-weighted average
#' @param plotKernel if TRUE, plots the kernel
#' @export
#' @examples
#' s = spectrogram(soundgen(), samplingRate = 16000, windowLength = 10,
#' output = 'original', plot = FALSE)
#' s = log(s + .001)
#' # image(s)
#' s1 = gaussianSmooth2D(s, kernelSize = 5, plotKernel = TRUE)
#' # image(s1)
#'
#' \dontrun{
#' # more smoothing in time than in frequency
#' s2 = gaussianSmooth2D(s, kernelSize = c(5, 15))
#' image(s2)
#'
#' # vice versa - more smoothing in frequency
#' s3 = gaussianSmooth2D(s, kernelSize = c(25, 3))
#' image(s3)
#'
#' # sharpen the image by deconvolution with the kernel
#' s4 = gaussianSmooth2D(s1, kernelSize = 5, action = 'unblur')
#' image(s4)
#'
#' s5 = gaussianSmooth2D(s, kernelSize = c(15, 1), action = 'unblur')
#' image(s5)
#' }
gaussianSmooth2D = function(m,
kernelSize = 5,
kernelSD = .5,
action = c('blur', 'unblur')[1],
plotKernel = FALSE) {
if (length(kernelSize) == 1) kernelSize = c(kernelSize, kernelSize)
nr = nrow(m)
nc = ncol(m)
if (max(kernelSize) < 2) return(m)
if (kernelSize[1] >= (nr / 2)) {
kernelSize[1] = ceiling(nr / 2) - 1
}
if (kernelSize[2] >= (nc / 2)) {
kernelSize[2] = ceiling(nc / 2) - 1
}
if (kernelSize[1] %% 2 == 0) kernelSize[1] = kernelSize[1] + 1 # make uneven
if (kernelSize[2] %% 2 == 0) kernelSize[2] = kernelSize[2] + 1 # make uneven
if (max(kernelSize) < 2) return(m)
# set up 2D Gaussian filter
kernel = getCheckerboardKernel(
size = kernelSize,
kernelSD = kernelSD,
checker = FALSE,
plot = plotKernel)
kernel = kernel / sum(kernel) # convert to pdf
## pad matrix with size / 2 zeros, so that we can correlate it with the
# kernel starting from the very edge
m_padded = matrix(0,
nrow = nr + kernelSize[1],
ncol = nc + kernelSize[2])
# lower left corner in the padded matrix where we'll paste the original m
idx_row = round((kernelSize[1] + 1) / 2)
idx_col = round((kernelSize[2] + 1) / 2)
# paste original. Now we have a padded matrix
m_padded[idx_row:(idx_row + nr - 1), idx_col:(idx_col + nc - 1)] = m
# kernel smoothing / deconvolution
out = matrix(NA, nrow = nr, ncol = nc)
if (action == 'blur') {
for (i in 1:nr) {
for (j in 1:nc) {
out[i, j] = sum(m_padded[i : (i + 2 * idx_row - 2), j : (j + 2 * idx_col - 2)] * kernel)
}
}
} else if (action == 'unblur') {
for (i in 1:nr) {
for (j in 1:nc) {
out[i, j] = m[i, j] - sum(m_padded[i : (i + 2 * idx_row - 2), j : (j + 2 * idx_col - 2)] * kernel)
}
}
}
if (!is.null(rownames(m))) rownames(out) = rownames(m)
if (!is.null(colnames(m))) colnames(out) = colnames(m)
return(out)
}
#' Proportion of total
#'
#' Internal soundgen function.
#'
#' Calculates the values in the input distribution that contain particular
#' proportions of the sum of all values in the input distribution.
#' @param x numeric vector of non-negative real numbers
#' @param quantiles quantiles of the cumulative distribution
#' @keywords internal
#' @examples
#' x = rnorm(100)
#' x = x - min(x) # must be non-negative
#' hist(x)
#' v = soundgen:::pDistr(x, quantiles = c(.5, .8, .9))
#' sum(x[x > v['0.5']]) / sum(x)
#' sum(x[x > v['0.9']]) / sum(x)
pDistr = function(x, quantiles) {
a1 = rev(sort(x)) # plot(a1, type = 'l')
a2 = cumsum(a1) # plot(a2, type = 'l')
total = sum(x)
out = rep(NA, length(quantiles))
names(out) = quantiles
for (q in 1:length(quantiles)) {
out[q] = a1[which(a2 > total * quantiles[q])[1]]
}
return(out)
}
#' Log-warp matrix
#'
#' Internal soundgen function.
#'
#' Log-warps a matrix, as if log-transforming plot axes.
#'
#' @param m a matrix of numeric values of any dimensions (not necessarily
#' square)
#' @param base the base of logarithm
#' @keywords internal
#' @examples
#' m = matrix(1:90, nrow = 10)
#' colnames(m) = 1:9
#' soundgen:::logMatrix(m, base = 2)
#' soundgen:::logMatrix(m, base = 10)
#'
#' soundgen:::logMatrix(m = matrix(1:9, nrow = 1), base = 2)
#'
#' \dontrun{
#' s = spectrogram(soundgen(), 16000, output = 'original')
#' image(log(t(soundgen:::logMatrix(s, base = 2))))
#' }
logMatrix = function(m, base = 2) {
# the key is to make a sequence of time locations within each row/column for
# interpolation: (1:nrow(m)) ^ base (followed by normalization, so this index
# will range from 1 to nrow(m) / ncol(m))
if (ncol(m) > 1) {
idx_row = (1:ncol(m)) ^ base - 1
idx_row = idx_row / max(idx_row)
idx_row = idx_row * (ncol(m) - 1) + 1
# interpolate rows at these time points
m1 = t(apply(m, 1, function(x) approx(x, xout = idx_row)$y))
} else {
m1 = m
}
# same for columns
if (nrow(m) > 1) {
idx_col = (1:nrow(m)) ^ base - 1
idx_col = idx_col / max(idx_col)
idx_col = idx_col * (nrow(m) - 1) + 1
# interpolate columns at these time points
m2 = apply(m1, 2, function(x) approx(x, xout = idx_col)$y)
} else {
m2 = m1
}
# interpolate row and column names
# (assuming numeric values, as when called from modulationSpectrum())
if (!is.null(colnames(m)) & ncol(m) > 1) {
colnames(m2) = approx(as.numeric(colnames(m)), xout = idx_row)$y
}
if (!is.null(rownames(m)) & nrow(m) > 1) {
rownames(m2) = approx(as.numeric(rownames(m)), xout = idx_col)$y
}
return(m2)
}
#' Switch color theme
#'
#' Internal soundgen function
#' @param colorTheme string like 'bw', 'seewave', or function name
#' @keywords internal
switchColorTheme = function(colorTheme) {
if (is.null(colorTheme)) return(NULL)
if (colorTheme == 'bw') {
color.palette = function(x) gray(seq(from = 1, to = 0, length = x))
} else if (colorTheme == 'seewave') {
color.palette = seewave::spectro.colors
} else {
colFun = match.fun(colorTheme)
color.palette = function(x) rev(colFun(x))
}
return(color.palette)
}
#' Parabolic peak interpolation
#'
#' Internal soundgen function
#'
#' Takes a spectral peak and two adjacent points, fits a parabola through them,
#' and thus estimates true peak location relative to the discrete peak. See
#' https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html
#' @param threePoints the amplitudes of three adjacent points (beta is the
#' peak), ideally spectrum on a dB scale, obtained with a Gaussian window
#' @param plot if TRUE, plot the points and the fit parabola
#' @return Returns a list: $p = the correction coefficient in bins (idx_beta + p
#' gives the true peak), $ampl_p = the amplitude of the true peak
#' @keywords internal
#' @examples
#' soundgen:::parabPeakInterpol(c(-1, 0, -4), plot = TRUE)
parabPeakInterpol = function(threePoints, plot = FALSE) {
if (any(is.na(threePoints)) | length(threePoints) < 3) {
stop(paste('invalid input:', threePoints))
}
if (threePoints[2] < threePoints[1] |
threePoints[2] < threePoints[3]) {
# the central point is not a peak
return(list(p = 0, ampl_p = 0))
}
alpha = threePoints[1]
beta = threePoints[2]
gamma = threePoints[3]
a = (alpha - 2 * beta + gamma) / 2
p = (alpha - gamma) / 4 / a
ampl_p = beta - (alpha - gamma) * p / 4
if (plot) {
b = beta - a * p ^ 2
x = seq(-2, 2, length.out = 100)
parab = a * (x - p) ^ 2 + b
plot(-1:1, c(alpha, beta, gamma), type = 'p',
xlim = c(-2, 2), ylim = c(gamma, ampl_p),
xlab = 'Bin', ylab = 'Ampl')
points(x, parab, type = 'l', col = 'blue')
}
return(list(p = p, ampl_p = ampl_p))
}
#' Identify and play
#'
#' Internal soundgen function
#'
#' A wrapper around \code{identify()} intended to play the sound corresponding
#' to a clicked point.
#' @param x,y plot coordinates
#' @param data dataframe from which x & y are taken, also containing a column
#' called "file"
#' @param audioFolder path to audio files
#' @param to play only the first ... s
#' @param plot if TRUE, plots the index of clicked points
#' @param pch symbol for marking clicked points
#' @param ... other arguments passed to \code{identify()}
#' @keywords internal
#' @examples
#' \dontrun{
#' msf = modulationSpectrum('~/Downloads/temp', plot = FALSE)
#'
#' # Method 1: provide path to folder, leave data = NULL
#' plot(msf$summary$amFreq_median, msf$summary$amDep_median)
#' identifyPch(msf$summary$amFreq_median, msf$summary$amDep_median,
#' audioFolder = '~/Downloads/temp',
#' to = 2,
#' plot = TRUE,
#' pch = 19)
#'
#' # Method 2:
#' plot(msf$summary$amFreq_median, msf$summary$amDep_median)
#' x = msf$summary$amFreq_median
#' y = msf$summary$amDep_median
#' identifyPch(x, y, data = msf$summary,
#' audioFolder = '~/Downloads/temp',
#' to = 2,
#' plot = FALSE,
#' pch = 8)
#' }
identifyAndPlay = function(
x,
y = NULL,
data = NULL,
audioFolder,
to = 5,
plot = FALSE,
pch = 19,
...) {
n = length(x)
if (is.null(data)) {
data = data.frame(
file = list.files(audioFolder, full.names = TRUE)
)
} else {
data = data.frame(
file = paste0(audioFolder, '/', data$file)
)
}
xy = xy.coords(x, y)
x = xy$x
y = xy$y
sel = rep(FALSE, n)
answers = numeric(0)
while(sum(sel) < n) {
ans = identify(x[!sel], y[!sel], labels = which(!sel),
n = 1, plot = plot, ...)
if(!length(ans)) break
ans = which(!sel)[ans]
answers = c(answers, ans)
points(x[ans], y[ans], pch = pch)
## play the selected point
try(playme(data$file[ans], to = to))
}
## return indices of selected points
return(data$file[answers])
}
#' Warp matrix
#'
#' Internal soundgen function
#'
#' Warps or scales each column of a matrix (normally a spectrogram).
#' @keywords internal
#' @param m matrix (rows = frequency bins, columns = time)
#' @param scaleFactor 1 = no change, >1 = raise formants
#' @param interpol interpolation method
#' @examples
#' a = matrix(1:12, nrow = 4)
#' a
#' soundgen:::warpMatrix(a, 1.5, 'approx')
#' soundgen:::warpMatrix(a, 1/1.5, 'spline')
warpMatrix = function(m, scaleFactor, interpol = c('approx', 'spline')[1]) {
scaleFactor = getSmoothContour(scaleFactor, len = ncol(m))
n1 = nrow(m)
m_warped = m
for (i in 1:ncol(m)) {
if (scaleFactor[i] > 1) {
# "stretch" the vector (eg spectrum of a frame)
n2 = round(n1 / scaleFactor[i])
m_warped[, i] = do.call(interpol, list(x = m[1:n2, i], n = n1))$y
} else if (scaleFactor[i] < 1) {
# "shrink" the vector and pad it with the last obs to the original length
n2 = round(n1 * scaleFactor[i])
padding = rep(m[n1, i], n1 - n2) # or 0
m_warped[, i] = c(
do.call(interpol, list(x = m_warped[, i],
xout = seq(1, n1, length.out = n2),
n = n1))$y,
padding
)
}
# plot(m_warped[, i], type = 'l')
# plot(abs(m[, i]), type = 'l', xlim = c(1, max(n1, n2)))
# points(m_warped[, i], type = 'l', col = 'blue')
}
return(m_warped)
}
#' Principal argument
#'
#' Internal soundgen function.
#'
#' Recalculates the phase of complex numbers to be within the interval from -pi to pi.
#' @param x real number representing phase angle
#' @keywords internal
princarg = function(x) (x + pi) %% (2 * pi) - pi
#' Split vector into chunks
#'
#' Internal soundgen function.
#'
#' Takes a numeric vector x and splits it into n chunks. This is the fastest
#' splitting algorithm from
#' https://stackoverflow.com/questions/3318333/split-a-vector-into-chunks
#' @param x numeric vector to split
#' @param n number of chunks
#' @return Returns a list of length \code{n} containing the chunks
#' @keywords internal
#' @examples
#' # prepare chunks of iterator to run in parallel on several cores
#' chunks = soundgen:::splitIntoChunks(1:21, 4)
#' chunks
splitIntoChunks = function(x, n) {
len = length(x)
len_chunk = ceiling(len / n)
mapply(function(a, b) (x[a:b]),
seq.int(from = 1, to = len, by = len_chunk),
pmin(seq.int(from = 1, to = len, by = len_chunk) + (len_chunk - 1), len),
SIMPLIFY = FALSE)
}
#' rbind_fill
#'
#' Internal soundgen function
#'
#' Fills missing columns with NAs, then rbinds - handy in case one has extra
#' columns. Used in formant_app(), pitch_app()
#' @param df1,df2 two dataframes with partly matching columns
#' @keywords internal
rbind_fill = function(df1, df2) {
if (!is.list(df1) || nrow(df1) == 0) return(df2)
if (!is.list(df2) || nrow(df2) == 0) return(df1)
df1[setdiff(names(df2), names(df1))] = NA
df2[setdiff(names(df1), names(df2))] = NA
return(rbind(df1, df2))
}
#' Pseudolog
#'
#' Internal soundgen function
#'
#' From function pseudo_log_trans() in the "scales" package.
#'
#' @seealso \code{\link{pseudoLog_undo}}
#'
#' @param x numeric vector to transform
#' @param sigma scaling factor for the linear part
#' @param base approximate logarithm base used
#' @keywords internal
#' @examples
#' a = -30:30
#' plot(a, soundgen:::pseudoLog(a, sigma = 1, base = 2))
#' plot(a, soundgen:::pseudoLog(a, sigma = 5, base = 2))
#' plot(a, soundgen:::pseudoLog(a, sigma = .1, base = 2))
pseudoLog = function(x, sigma, base) {
asinh(x/(2 * sigma))/log(base)
}
#' Undo pseudolog
#'
#' Internal soundgen function
#'
#' From function pseudo_log_trans() in the "scales" package.
#'
#' @seealso \code{\link{pseudoLog}}
#'
#' @param x numeric vector to transform
#' @param sigma scaling factor for the linear part
#' @param base approximate logarithm base used
#' @keywords internal
pseudoLog_undo = function(x, sigma, base) {
2 * sigma * sinh(x * log(base))
}
#' Find the elbow of a screeplot or similar
#'
#' Adapted from
#' https://stackoverflow.com/questions/2018178/finding-the-best-trade-off-point-on-a-curve
#' Algorithm: draw a straight line between the two endpoints and find the point
#' furthest from this line.
#'
#' @param d dataframe containing x and y coordinates of the points
#' @keywords internal
#' @examples
#' y = c(10, 11, 8, 4, 2, 1.5, 1, 0.7, .5, .4, .3)
#' soundgen:::findElbow(data.frame(x = 1:length(y), y = y), plot = TRUE)
findElbow = function(d, plot = FALSE) {
d = na.omit(d)
n = nrow(d)
# draw a straight line between the endpoints
endpoints = d[c(1, n), ]
fit = lm(endpoints$y ~ endpoints$x)
# calculate distances from points to line
# https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line
coefs = coef(fit)
distances = abs(coefs[2] * d$x - d$y + coefs[1]) / sqrt(coefs[2]^2 + 1)
out = which.max(distances)
if (plot) {
plot(d, type = 'b')
abline(coefs, col = 'blue', lty = 2)
# segment to the closest point on the line (again from wiki)
a = -coefs[2]; b = 1; c = -coefs[1]
x0 = d$x[out]; y0 = d$y[out]
x1 = (b * (b * x0 - a * y0) - a * c) / (a^2 + b^2)
y1 = (a * (-b * x0 + a * y0) - b * c) / (a^2 + b^2)
segments(x0, y0, x1, y1, lty = 3)
# NB: the angle only looks straight if the plot is square!
}
return(out)
}
#' Logistic
#' @param x numeric vector
#' @keywords internal
logistic = function(x) 1 / (1 + exp(-x))
#' Logit
#' @param x numeric vector
#' @keywords internal
logit = function(x) log(x / (1 - x))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.