Nothing
R/modulationSpectrum.R
In soundgen: Sound Synthesis and Acoustic Analysis
Defines functions
logWarpMS
getAM
averageMatrices
getRough
modulationSpectrumFragment
plotMS
.modulationSpectrum
modulationSpectrum
Documented in
averageMatrices
getAM
getRough
logWarpMS
modulationSpectrum
.modulationSpectrum
modulationSpectrumFragment
plotMS
### MODULATION SPECTRUM
#' Modulation spectrum
#'
#' Produces a modulation spectrum of waveform(s) or audio file(s), with temporal
#' modulation along the X axis (Hz) and spectral modulation (1/KHz) along the Y
#' axis. A good visual analogy is decomposing the spectrogram into a sum of
#' ripples of various frequencies and directions. Roughness is calculated as the
#' proportion of energy / amplitude of the modulation spectrum within
#' \code{roughRange} of temporal modulation frequencies. The frequency of
#' amplitude modulation (amMsFreq, Hz) is calculated as the highest peak in the
#' smoothed AM function, and its purity (amMsPurity, dB) as the ratio of this
#' peak to the median AM over \code{amRange}. For relatively short and steady
#' sounds, set \code{amRes = NULL} and analyze the entire sound. For longer
#' sounds and when roughness or AM vary over time, set \code{amRes} to get
#' multiple measurements over time (see examples).
#'
#' Algorithm: prepare a spectrogram, take its logarithm (if \code{logSpec =
#' TRUE}), center, perform a 2D Fourier transform (see also
#' spectral::spec.fft()), take the upper half of the resulting symmetric matrix,
#' and raise it to \code{power}. The result is returned as \code{$original}. For
#' plotting purposes, the modulation matrix can be smoothed with Gaussian blur
#' (see \code{\link{gaussianSmooth2D}}) and log-warped (if \code{logWarp} is a
#' positive number). This processed modulation spectrum is returned as
#' \code{$processed}. If the audio is long enough, multiple windows are
#' analyzed, resulting in a vector of roughness values. For multiple inputs,
#' such as a list of waveforms or path to a folder with audio files, the
#' ensemble of modulation spectra can be interpolated to the same spectral and
#' temporal resolution and averaged (if \code{averageMS}).
#'
#' @seealso \code{\link{plotMS}} \code{\link{spectrogram}} \code{\link{analyze}}
#'
#' @references \itemize{
#' \item Singh, N. C., & Theunissen, F. E. (2003). Modulation spectra of
#' natural sounds and ethological theories of auditory processing. The Journal
#' of the Acoustical Society of America, 114(6), 3394-3411.
#' }
#' @inheritParams spectrogram
#' @inheritParams analyze
#' @param roughRange the range of temporal modulation frequencies that
#' constitute the "roughness" zone, Hz
#' @param amRange the range of temporal modulation frequencies that we are
#' interested in as "amplitude modulation" (AM), Hz
#' @param amRes target resolution of amplitude modulation, Hz. If \code{NULL},
#' the entire sound is analyzed at once, resulting in a single roughness value
#' (unless it is longer than \code{maxDur}, in which case it is analyzed in
#' chunks \code{maxDur} s long). If \code{amRes} is set, roughness is
#' calculated for windows \code{~1000/amRes} ms long (but at least 3 STFT
#' frames). \code{amRes} also affects the amount of smoothing when calculating
#' \code{amMsFreq} and \code{amMsPurity}
#' @param maxDur sounds longer than \code{maxDur} s are split into fragments,
#' and the modulation spectra of all fragments are averaged
#' @param logSpec if TRUE, the spectrogram is log-transformed prior to taking 2D
#' FFT
#' @param power raise modulation spectrum to this power (eg power = 2 for ^2, or
#' "power spectrum")
#' @param normalize if TRUE, the modulation spectrum of each analyzed fragment
#' \code{maxDur} in duration is separately normalized to have max = 1
#' @param returnMS if FALSE, only roughness is returned (much faster)
#' @param returnComplex if TRUE, returns a complex modulation spectrum (without
#' normalization and warping)
#' @param averageMS if TRUE, the modulation spectra of all inputs are averaged
#' into a single output; if FALSE, a separate MS is returned for each input
#' @param plot if TRUE, plots the modulation spectrum of each sound
#' @param savePlots if a valid path is specified, a plot is saved in this folder
#' (defaults to NA)
#' @param logWarpX,logWarpY numeric vector of length 2: c(sigma, base) of
#' pseudolog-warping the modulation spectrum, as in function
#' pseudo_log_trans() from the "scales" package
#' @param quantiles labeled contour values, \% (e.g., "50" marks regions that
#' contain 50\% of the sum total of the entire modulation spectrum)
#' @param kernelSize the size of Gaussian kernel used for smoothing (1 = no
#' smoothing)
#' @param kernelSD the SD of Gaussian kernel used for smoothing, relative to its
#' size
#' @param xlab,ylab,main,xlim,ylim graphical parameters
#' @param width,height,units,res parameters passed to
#' \code{\link[grDevices]{png}} if the plot is saved
#' @param ... other graphical parameters passed on to \code{filled.contour.mod}
#' and \code{\link[graphics]{contour}} (see \code{\link{spectrogram}})
#' @return Returns a list with the following components:
#' \itemize{
#' \item \code{$original} modulation spectrum prior to blurring and log-warping,
#' but after squaring if \code{power = TRUE}, a matrix of nonnegative values.
#' Rownames are spectral modulation frequencies (cycles/KHz), and colnames are
#' temporal modulation frequencies (Hz).
#' \item \code{$processed} modulation spectrum after blurring and log-warping
#' \item \code{$complex} untransformed complex modulation spectrum (returned
#' only if returnComplex = TRUE)
#' \item \code{$roughness} proportion of energy / amplitude of the modulation
#' spectrum within \code{roughRange} of temporal modulation frequencies, \% - a
#' vector if amRes is numeric and the sound is long enough, a single number
#' otherwise
#' \item \code{$amMsFreq} frequency of the highest peak, within \code{amRange}, of
#' the folded AM function (average AM across all FM bins for both negative and
#' positive AM frequencies), where a peak is a local maximum over \code{amRes}
#' Hz. Like \code{roughness}, \code{amMsFreq} and \code{amMsPurity} can be single
#' numbers or vectors, depending on whether the sound is analyzed as a whole or
#' in chunks
#' \item \code{$amMsPurity} ratio of the peak at amMsFreq to the median AM over
#' \code{amRange}, dB
#' \item \code{$summary} dataframe with summaries of roughness, amMsFreq, and
#' amMsPurity
#' }
#' @export
#' @examples
#' # White noise
#' ms = modulationSpectrum(runif(16000), samplingRate = 16000,
#' logSpec = FALSE, power = TRUE,
#' amRes = NULL) # analyze the entire sound, giving a single roughness value
#' str(ms)
#'
#' # Harmonic sound
#' s = soundgen(amFreq = 25, amDep = 50)
#' ms = modulationSpectrum(s, samplingRate = 16000, amRes = NULL)
#' ms[c('roughness', 'amMsFreq', 'amMsPurity')] # a single value for each
#' ms1 = modulationSpectrum(s, samplingRate = 16000, amRes = 5)
#' ms1[c('roughness', 'amMsFreq', 'amMsPurity')]
#' # measured over time (low values of amRes mean more precision, so we analyze
#' # longer segments and get fewer values per sound)
#'
#' # Embellish
#' ms = modulationSpectrum(s, samplingRate = 16000,
#' xlab = 'Temporal modulation, Hz', ylab = 'Spectral modulation, 1/KHz',
#' colorTheme = 'heat.colors', main = 'Modulation spectrum', lty = 3)
#'
#' \dontrun{
#' # A long sound with varying AM and a bit of chaos at the end
#' s_long = soundgen(sylLen = 3500, pitch = c(250, 320, 280),
#' amFreq = c(30, 55), amDep = c(20, 60, 40),
#' jitterDep = c(0, 0, 2))
#' playme(s_long)
#' ms = modulationSpectrum(s_long, 16000)
#' # plot AM over time
#' plot(x = seq(1, 1500, length.out = length(ms$amMsFreq)), y = ms$amMsFreq,
#' cex = 10^(ms$amMsPurity/20) * 10, xlab = 'Time, ms', ylab = 'AM frequency, Hz')
#' # plot roughness over time
#' spectrogram(s_long, 16000, ylim = c(0, 4),
#' extraContour = list(ms$roughness / max(ms$roughness) * 4000, col = 'blue'))
#'
#' # As with spectrograms, there is a tradeoff in time-frequency resolution
#' s = soundgen(pitch = 500, amFreq = 50, amDep = 100,
#' samplingRate = 44100, plot = TRUE)
#' # playme(s, samplingRate = 44100)
#' ms = modulationSpectrum(s, samplingRate = 44100,
#' windowLength = 50, step = 50, amRes = NULL) # poor temporal resolution
#' ms = modulationSpectrum(s, samplingRate = 44100,
#' windowLength = 5, step = 1, amRes = NULL) # poor frequency resolution
#' ms = modulationSpectrum(s, samplingRate = 44100,
#' windowLength = 15, step = 3, amRes = NULL) # a reasonable compromise
#'
#' # customize the plot
#' ms = modulationSpectrum(s, samplingRate = 44100,
#' windowLength = 15, overlap = 80, amRes = NULL,
#' kernelSize = 17, # more smoothing
#' xlim = c(-70, 70), ylim = c(0, 4), # zoom in on the central region
#' quantiles = c(.25, .5, .8), # customize contour lines
#' col = rev(rainbow(100)), # alternative palette
#' power = 2) # ^2
#' # Note the peaks at FM = 2/KHz (from "pitch = 500") and AM = 50 Hz (from
#' # "amFreq = 50")
#'
#' # Input can be a wav/mp3 file
#' ms = modulationSpectrum('~/Downloads/temp/16002_Faking_It_Large_clear.wav')
#'
#' # Input can be path to folder with audio files. Each file is processed
#' # separately, and the output can contain an MS per file...
#' ms1 = modulationSpectrum('~/Downloads/temp', kernelSize = 11,
#' plot = FALSE, averageMS = FALSE)
#' ms1$summary
#' names(ms1$original) # a separate MS per file
#' # ...or a single MS can be calculated:
#' ms2 = modulationSpectrum('~/Downloads/temp', kernelSize = 11,
#' plot = FALSE, averageMS = TRUE)
#' plotMS(ms2$original)
#' ms2$summary
#'
#' # Input can also be a list of waveforms (numeric vectors)
#' ss = vector('list', 10)
#' for (i in 1:length(ss)) {
#' ss[[i]] = soundgen(sylLen = runif(1, 100, 1000), temperature = .4,
#' pitch = runif(3, 400, 600))
#' }
#' # lapply(ss, playme)
#' # MS of the first sound
#' ms1 = modulationSpectrum(ss[[1]], samplingRate = 16000, scale = 1)
#' # average MS of all 10 sounds
#' ms2 = modulationSpectrum(ss, samplingRate = 16000, scale = 1, averageMS = TRUE, plot = FALSE)
#' plotMS(ms2$original)
#'
#' # A sound with ~3 syllables per second and only downsweeps in F0 contour
#' s = soundgen(nSyl = 8, sylLen = 200, pauseLen = 100, pitch = c(300, 200))
#' # playme(s)
#' ms = modulationSpectrum(s, samplingRate = 16000, maxDur = .5,
#' xlim = c(-25, 25), colorTheme = 'seewave',
#' power = 2)
#' # note the asymmetry b/c of downsweeps
#'
#' # "power = 2" returns squared modulation spectrum - note that this affects
#' # the roughness measure!
#' ms$roughness
#' # compare:
#' modulationSpectrum(s, samplingRate = 16000, maxDur = .5,
#' xlim = c(-25, 25), colorTheme = 'seewave',
#' power = 1)$roughness # much higher roughness
#'
#' # Plotting with or without log-warping the modulation spectrum:
#' ms = modulationSpectrum(soundgen(), samplingRate = 16000, plot = TRUE)
#' ms = modulationSpectrum(soundgen(), samplingRate = 16000,
#' logWarpX = c(2, 2), plot = TRUE)
#'
#' # logWarp and kernelSize have no effect on roughness
#' # because it is calculated before these transforms:
#' modulationSpectrum(s, samplingRate = 16000, logWarpX = c(1, 10))$roughness
#' modulationSpectrum(s, samplingRate = 16000, logWarpX = NA)$roughness
#' modulationSpectrum(s, samplingRate = 16000, kernelSize = 17)$roughness
#'
#' # Log-transform the spectrogram prior to 2D FFT (affects roughness):
#' modulationSpectrum(s, samplingRate = 16000, logSpec = FALSE)$roughness
#' modulationSpectrum(s, samplingRate = 16000, logSpec = TRUE)$roughness
#'
#' # Complex modulation spectrum with phase preserved
#' ms = modulationSpectrum(soundgen(), samplingRate = 16000,
#' returnComplex = TRUE)
#' plotMS(abs(ms$complex)) # note the symmetry
#' # compare:
#' plotMS(ms$original)
#' }
modulationSpectrum = function(
x,
samplingRate = NULL,
scale = NULL,
from = NULL,
to = NULL,
amRes = 5,
maxDur = 5,
logSpec = FALSE,
windowLength = 15,
step = NULL,
overlap = 80,
wn = 'hanning',
zp = 0,
power = 1,
normalize = TRUE,
roughRange = c(30, 150),
amRange = c(10, 200),
returnMS = TRUE,
returnComplex = FALSE,
summaryFun = c('mean', 'median', 'sd'),
averageMS = FALSE,
reportEvery = NULL,
cores = 1,
plot = TRUE,
savePlots = NULL,
logWarpX = NULL,
logWarpY = NULL,
quantiles = c(.5, .8, .9),
kernelSize = 5,
kernelSD = .5,
colorTheme = c('bw', 'seewave', 'heat.colors', '...')[1],
col = NULL,
main = NULL,
xlab = 'Hz',
ylab = '1/KHz',
xlim = NULL,
ylim = NULL,
width = 900,
height = 500,
units = 'px',
res = NA,
...
) {
## Prepare a list of arguments to pass to .modulationSpectrum()
myPars = c(as.list(environment()), list(...))
# exclude unnecessary args
myPars = myPars[!names(myPars) %in% c(
'x', 'samplingRate', 'scale', 'from', 'to', 'savePlots',
'reportEvery', 'cores', 'summaryFun', 'averageMS')]
# analyze
pa = processAudio(
x,
samplingRate = samplingRate,
scale = scale,
from = from,
to = to,
funToCall = '.modulationSpectrum',
myPars = myPars,
reportEvery = reportEvery,
cores = cores,
savePlots = savePlots
)
# htmlPlots
if (!is.null(pa$input$savePlots) && pa$input$n > 1) {
try(htmlPlots(pa$input, savePlots = savePlots, changesAudio = FALSE,
suffix = "MS", width = paste0(width, units)))
}
# prepare output
if (!is.null(summaryFun) && any(!is.na(summaryFun))) {
temp = vector('list', pa$input$n)
for (i in 1:pa$input$n) {
if (!pa$input$failed[i]) {
temp[[i]] = summarizeAnalyze(
data.frame(roughness = pa$result[[i]]$roughness,
amMsFreq = pa$result[[i]]$amMsFreq,
amMsPurity = pa$result[[i]]$amMsPurity),
summaryFun = summaryFun,
var_noSummary = NULL)
}
}
idx_failed = which(pa$input$failed)
idx_ok = which(!pa$input$failed)
if (length(idx_failed) > 0) {
idx_ok = which(!pa$input$failed)
if (length(idx_ok) > 0) {
filler = temp[[idx_ok[1]]] [1, ]
filler[1, ] = NA
} else {
stop('Failed to analyze any input')
}
for (i in idx_failed) temp[[i]] = filler
}
mysum_all = cbind(data.frame(file = pa$input$filenames_base),
do.call('rbind', temp))
} else {
mysum_all = NULL
}
for (i in idx_failed) pa$result[[i]] = list(
original = NULL, processed = NULL, complex = NULL,
roughness = NULL, amMsFreq = NULL, amMsPurity = NULL
)
original = processed = complex = roughness = amMsFreq = amMsPurity = NULL
# (otherwise note about no visible binding)
out_prep = c('original', 'processed', 'complex', 'roughness', 'amMsFreq', 'amMsPurity')
if (pa$input$n == 1) {
# unlist
for (op in out_prep) assign(noquote(op), pa$result[[1]] [[op]])
} else {
for (op in out_prep) assign(noquote(op), lapply(pa$result, function(x) x[[op]]))
}
# average MS across sounds
op_prep2 = c('original', 'processed') # roughness is not a matrix
if (returnComplex) {
op_prep2 = c(op_prep2, 'complex')
} else {
complex = NULL
}
if (averageMS & pa$input$n > 1) {
for (op in op_prep2) {
assign(
noquote(op),
averageMatrices(
get(op) [idx_ok],
rFun = 'max', # normally same samplingRate, but if not, upsample frequency resolution
cFun = 'median', # typical ncol (depends on sound dur)
reduceFun = '+') # internally divides by length, so actually becomes mean
)
}
}
invisible(list(original = original,
processed = processed,
complex = complex,
roughness = roughness,
amMsFreq = amMsFreq,
amMsPurity = amMsPurity,
summary = mysum_all))
}
#' Modulation spectrum per sound
#'
#' Internal soundgen function.
#' @inheritParams modulationSpectrum
#' @param audio a list returned by \code{readAudio}
#' @keywords internal
.modulationSpectrum = function(
audio,
amRes = 5,
maxDur = 5,
logSpec = FALSE,
windowLength = 15,
step = NULL,
overlap = 80,
wn = 'hanning',
zp = 0,
power = 1,
normalize = TRUE,
roughRange = c(30, 150),
amRange = c(10, 200),
returnMS = TRUE,
returnComplex = FALSE,
plot = TRUE,
savePlots = NULL,
logWarpX = NULL,
logWarpY = NULL,
quantiles = c(.5, .8, .9),
kernelSize = 5,
kernelSD = .5,
colorTheme = c('bw', 'seewave', 'heat.colors', '...')[1],
col = NULL,
main = NULL,
xlab = 'Hz',
ylab = '1/KHz',
xlim = NULL,
ylim = NULL,
width = 900,
height = 500,
units = 'px',
res = NA,
...
) {
# Re-set windowLength, step, and overlap so as to ensure that
# windowLength_points and step_points are not fractions
if (is.null(step)) step = windowLength * (1 - overlap / 100)
step_points = round(step / 1000 * audio$samplingRate)
step = step_points / audio$samplingRate * 1000
windowLength_points = round(windowLength / 1000 * audio$samplingRate)
windowLength = windowLength_points / audio$samplingRate * 1000
overlap = 100 * (1 - step_points / windowLength_points)
max_am = 1000 / step / 2
if (max_am < roughRange[1]) {
warning(paste(
'roughRange outside the analyzed range of temporal modulation frequencies;',
'increase overlap / decrease step to improve temporal resolution,',
'or else look for roughness in a lower range'))
}
if (is.numeric(amRes)) {
nFrames = max(3, ceiling(max_am / amRes * 2)) # min 3 spectrograms frames
} else {
nFrames = NULL
}
if (is.numeric(nFrames)) {
# split the input sound into chunks nFrames long
chunk_ms = windowLength + step * (nFrames - 1)
splitInto = max(1, ceiling(audio$duration * 1000 / chunk_ms))
if (chunk_ms > (audio$duration * 1000)) {
message(paste('The sound is too short to be analyzed with amRes =', amRes,
'Hz. Roughness is probably not measured correctly'))
chunk_ms = audio$duration
}
} else {
# split only those sounds that exceed maxDur
splitInto = max(1, ceiling(audio$duration / maxDur))
# so, for ex., if 2.1 times longer than maxDur, split into three
}
if (splitInto > 1) {
myseq = floor(seq(1, length(audio$sound), length.out = splitInto + 1))
myInput = vector('list', splitInto)
for (i in 1:splitInto) {
idx = myseq[i]:(myseq[i + 1])
myInput[[i]] = audio$sound[idx]
}
} else {
myInput = list(audio$sound)
}
# extract modulation spectrum per fragment
out = vector('list', splitInto)
if (returnComplex) {
out_complex = out
} else {
out_aggreg_complex = NULL
}
for (i in 1:splitInto) {
ms_i = modulationSpectrumFragment(
myInput[[i]],
samplingRate = audio$samplingRate,
windowLength = windowLength,
windowLength_points = windowLength_points,
step = step,
step_points = step_points,
wn = wn,
zp = zp,
logSpec = logSpec,
power = power,
normalize = normalize)
out[[i]] = ms_i$ms_half
if (returnComplex) {
out_complex[[i]] = ms_i$ms_complex
}
}
keep = which(unlist(lapply(out, function(x) !is.null(x))))
out = out[keep]
if (length(out) < 1) {
warning('The sound is too short or windowLength too long. Need at least 3 frames')
return(list('original' = NA,
'processed' = NA,
'complex' = NA,
'roughness' = NA,
'amMsFreq' = NA,
'amMsPurity' = NA))
}
# extract a measure of roughness
roughness = unlist(lapply(out, function(x)
getRough(x, roughRange, colNames = as.numeric(colnames(out[[1]])))))
# detect systematic amplitude modulation at the same frequency
am_list = lapply(out, function(x)
getAM(x, amRange = amRange, amRes = amRes))
amMsFreq = unlist(lapply(am_list, function(x) x$amMsFreq))
amMsPurity = unlist(lapply(am_list, function(x) x$amMsPurity))
# average modulation spectra across all sounds
if (!returnMS) {
result = list('original' = NULL,
'processed' = NULL,
'complex' = NULL,
'roughness' = roughness,
'amMsFreq' = amMsFreq,
'amMsPurity' = amMsPurity)
} else {
out_aggreg = averageMatrices(
out,
rFun = 'max', # normally same samplingRate, but if not, upsample frequency resolution
cFun = 'median', # typical ncol (depends on sound dur)
reduceFun = '+') # internally divides by length, so actually becomes mean
X = as.numeric(colnames(out_aggreg))
Y = as.numeric(rownames(out_aggreg))
# image(X, Y, t(log(out_aggreg)))
# prepare a separate summary of the complex ms
if (returnComplex) {
out_aggreg_complex = averageMatrices(
out_complex,
rFun = 'max',
cFun = 'median',
reduceFun = '+')
# image(t(log(abs(out_aggreg_complex))))
}
# smoothing / blurring
out_transf = gaussianSmooth2D(out_aggreg,
kernelSize = kernelSize,
kernelSD = kernelSD)
result = list('original' = out_aggreg,
'processed' = out_transf,
'complex' = out_aggreg_complex,
'roughness' = roughness,
'amMsFreq' = amMsFreq,
'amMsPurity' = amMsPurity)
} # end of if (returnMS)
# PLOTTING
if (is.character(audio$savePlots)) {
plot = TRUE
png(filename = paste0(audio$savePlots, audio$filename_noExt, "_MS.png"),
width = width, height = height, units = units, res = res)
}
if (plot) {
plotMS(ms = out_transf,
X = X,
Y = Y,
audio = audio,
colorTheme = colorTheme,
col = col,
logWarpX = logWarpX,
logWarpY = logWarpY,
xlab = xlab, ylab = ylab,
main = main,
xlim = xlim, ylim = ylim,
quantiles = quantiles,
...
)
if (is.character(audio$savePlots)) dev.off()
}
invisible(result)
}
#' Plot modulation spectrum
#'
#' Plots a single modulation spectrum returned by
#' \code{\link{modulationSpectrum}}. The result is the same as the plot produced
#' by \code{\link{modulationSpectrum}}, but calling \code{plotMS} is handy for
#' processed modulation spectra - for instance, for plotting the difference
#' between the modulation spectra of two sounds or groups of sounds.
#'
#' @inheritParams modulationSpectrum
#' @param ms modulation spectrum - a matrix with temporal modulation in columns
#' and spectral modulation in rows, as returned by
#' \code{\link{modulationSpectrum}}
#' @param X,Y rownames and colnames of \code{ms}, respectively
#' @param audio (internal) a list of audio attributes
#' @export
#' @examples
#' ms1 = modulationSpectrum(runif(4000), samplingRate = 16000, plot = TRUE)
#' plotMS(ms1$processed) # identical to above
#'
#' # compare two modulation spectra
#' ms2 = modulationSpectrum(soundgen(sylLen = 100, addSilence = 0),
#' samplingRate = 16000)
#' # ensure the two matrices have the same dimensions
#' ms2_resized = soundgen:::interpolMatrix(ms2$original,
#' nr = nrow(ms1$original), nc = ncol(ms1$original))
#' # plot the difference
#' plotMS(log(ms1$original / ms2_resized), quantile = NULL,
#' col = colorRampPalette(c('blue', 'yellow')) (50))
plotMS = function(
ms,
X = NULL,
Y = NULL,
quantiles = c(.5, .8, .9),
colorTheme = c('bw', 'seewave', 'heat.colors', '...')[1],
col = NULL,
logWarpX = NULL,
logWarpY = NULL,
main = NULL,
xlab = 'Hz',
ylab = '1/KHz',
xlim = NULL,
ylim = NULL,
audio = NULL,
...
) {
if (!is.null(col)) colorTheme = NULL
if (!is.null(colorTheme)) {
color.palette = switchColorTheme(colorTheme)
} else {
color.palette = NULL
}
if (is.null(X)) X = as.numeric(colnames(ms))
if (is.null(Y)) Y = as.numeric(rownames(ms))
if (is.null(xlim)) xlim = c(X[1], -X[1])
if (is.null(ylim)) ylim = range(Y)
if (is.null(main) & !is.null(audio$filename_noExt)) {
if (audio$filename_noExt == 'sound') {
main = ''
} else {
main = audio$filename_noExt
}
}
# plot with filled.contour.mod
ms = zeroOne(ms)
X1 = X
Y1 = Y
if (is.numeric(logWarpX) | is.numeric(logWarpY)) {
if (is.numeric(logWarpX)) {
X1 = pseudoLog(X, sigma = logWarpX[1], base = logWarpX[2])
# lab_x = pretty(X, n = 9, min.n = 7)
# at_x = pseudoLog(lab_x, sigma = logWarpX[1], base = logWarpX[2])
at_x = pretty(X1, n = 7)
lab_x = round(pseudoLog_undo(at_x, sigma = logWarpX[1], base = logWarpX[2]))
} else {
lab_x = at_x = pretty(X)
}
if (is.numeric(logWarpY)) {
Y1 = pseudoLog(Y, sigma = logWarpY[1], base = logWarpY[2])
lab_y = pretty(Y)
at_y = pseudoLog(lab_y, sigma = logWarpY[2], base = logWarpY[2])
} else {
lab_y = at_y = pretty(Y)
}
filled.contour.mod(X1, Y1, t(ms),
levels = seq(0, 1, length = 100),
color.palette = color.palette,
col = col,
xlab = xlab, ylab = ylab,
bty = 'n',
main = main, xaxt = 'n', yaxt = 'n',
# xlim = xlim, ylim = ylim,
...)
# add nicely labeled axes
axis(1, at = at_x, labels = lab_x, ...) # xpd = TRUE to extend beyond the plot
axis(2, at = at_y, labels = lab_y, ...)
lbl_Hz_pos = pretty(Y1)
lbls_Hz = round(1000 / lbl_Hz_pos)
axis(4, at = lbl_Hz_pos, labels = lbls_Hz, ...)
} else {
filled.contour.mod(X1, Y1, t(ms),
levels = seq(0, 1, length = 100),
color.palette = color.palette,
col = col,
xlab = xlab, ylab = ylab,
bty = 'n',
main = main,
# xlim = xlim, ylim = ylim,
...)
lbls = round(1000 / pretty(Y1))
lbl_pos = 1000 / lbls
axis(4, at = lbl_pos, labels = lbls, ...)
}
abline(v = 0, lty = 3)
# add contours
if (is.numeric(quantiles)) {
# qntls = quantile(out_aggreg, probs = quantiles) # could try HDI instead
qntls = pDistr(as.numeric(ms), quantiles = quantiles)
par(new = TRUE)
contour(x = X1, y = Y1, z = t(ms),
levels = qntls, labels = quantiles * 100,
xaxs = 'i', yaxs = 'i',
axes = FALSE, frame.plot = FALSE,
# xlim = xlim, ylim = ylim,
...)
par(new = FALSE)
}
}
#' Modulation spectrum per fragment
#'
#' Internal soundgen function.
#' @inheritParams modulationSpectrum
#' @param sound numeric vector
#' @keywords internal
#' @examples
#' s = soundgen(amFreq = 25, amDep = 100)
#' ms = soundgen:::modulationSpectrumFragment(s, 16000,
#' windowLength = 50, windowLength_points = .05 * 16000,
#' step = 5, step_points = .005 * 16000)
#' image(as.numeric(colnames(ms$ms_half)), as.numeric(rownames(ms$ms_half)),
#' t(log(ms$ms_half)))
modulationSpectrumFragment = function(sound,
samplingRate,
windowLength,
windowLength_points,
step,
step_points,
wn = 'hanning',
zp = 0,
logSpec = FALSE,
power = 1,
normalize = TRUE) {
# Calling stdft is ~80 times faster than going through spectrogram (!)
step_seq = seq(1, length(sound) + 1 - windowLength_points, step_points)
# print(length(step_seq))
if (length(step_seq) < 3) return(NULL)
s1 = seewave::stdft(
wave = as.matrix(sound),
wn = wn,
wl = windowLength_points,
f = samplingRate,
zp = zp,
step = step_seq,
scale = FALSE,
norm = FALSE,
complex = FALSE
)
# image(t(log(s1))); s1[1:3, 1:3]; dim(s1); range(s1)
# log-transform amplitudes
if (logSpec) {
positives = which(s1 > 0)
nonpositives = which(s1 <= 0)
s1[positives] = log(s1[positives])
if (length(positives) > 0 & length(nonpositives) > 0) {
s1[nonpositives] = min(s1[positives])
}
s1 = s1 - min(s1) + 1e-16 # positive
# image(t(log(s1)))
}
# 2D fft
ms_complex = specToMS(s1, windowLength = windowLength, step = step)
ms = abs(ms_complex)
# image(as.numeric(colnames(ms)), as.numeric(rownames(ms)), t(log(ms)))
symAxis = floor(nrow(ms) / 2) + 1
# ms[(symAxis - 2) : (symAxis + 2), 1:2]
ms_half = ms[symAxis:nrow(ms),, drop = FALSE] # take only the upper half (always even)
# power
if (is.numeric(power) && power != 1) ms_half = ms_half ^ power
# normalize
if (normalize && any(s1 != 0)) {
ms_half = ms_half - min(ms_half)
ms_half = ms_half / max(ms_half)
}
# image(x = as.numeric(colnames(ms_half)), z = t(log(ms_half)))
return(list(
ms_half = ms_half,
ms_complex = ms_complex
))
}
#' Get roughness
#'
#' Internal soundgen function
#'
#' Helper function for calculating roughness - the proportion of energy /
#' amplitude in the roughness range
#' @param m numeric matrix of non-negative values with colnames giving temporal
#' modulation frequency
#' @inheritParams modulationSpectrum
#' @param colNames numeric vector with frequencies corresponding to columns
#' @return Returns roughness in percent.
#' @keywords internal
#' @examples
#' m = modulationSpectrum(soundgen(), samplingRate = 16000)$original
#' soundgen:::getRough(m, roughRange = c(30, Inf))
getRough = function(m, roughRange, colNames = NULL) {
if (is.null(colNames)) colNames = abs(as.numeric(colnames(m)))
rough_cols = which(colNames > roughRange[1] &
colNames < roughRange[2])
if (length(rough_cols) > 0) {
roughness = sum(m[, rough_cols]) / sum(m) * 100
} else {
roughness = 0
}
# # alternatively / in addition, could try gaussian filter when calculating roughness
# ampl_filter = dnorm(colNames, mean = roughMean, sd = roughSD)
# ampl_filter = ampl_filter / sum(ampl_filter) # to pdf
# # plot(ampl_filter, type = 'l')
# roughness = sum(colSums(m) * ampl_filter) / sum(m) * 100
return(roughness)
}
#' Average matrices
#'
#' Internal soundgen function.
#'
#' Takes a list of matrices (normally modulation spectra), interpolates them to
#' have the same size, and then reduces them (eg takes the average).
#' @param mat_list a list of matrices to aggregate (eg spectrograms or
#' modulation spectra)
#' @param rFun, cFun functions used to determine the number of rows and columns
#' in the result
#' @param reduceFun function used to aggregate
#' @keywords internal
#' @examples
#' mat_list = list(
#' matrix(1:30, nrow = 5),
#' matrix(80:17, nrow = 8)
#' )
#' soundgen:::averageMatrices(mat_list)
#' soundgen:::averageMatrices(mat_list, cFun = 'max', reduceFun = '*')
averageMatrices = function(mat_list,
rFun = 'max',
cFun = 'median',
reduceFun = '+') {
# normally same samplingRate, but in case not, upsample frequency resolution
nr = round(do.call(rFun, list(unlist(lapply(mat_list, nrow)))))
# take typical ncol (depends on sound dur)
nc = round(do.call(cFun, list(unlist(lapply(mat_list, ncol)))))
mat_list_sameDim = lapply(mat_list, function(x) interpolMatrix(x, nr = nr, nc = nc))
agg = Reduce(reduceFun, mat_list_sameDim) / length(mat_list)
return(agg)
}
#' Get amplitude modulation
#'
#' Internal soundgen function
#'
#' Helper function for calculating amplitude modulation based on the modulation
#' spectrum. Algorithm: averages AM across all FM bins in the positive half of
#' the modulation spectrum and looks for a peak in the specified AM frequency
#' range.
#' @param m numeric matrix of non-negative values with colnames giving temporal
#' modulation frequency
#' @inheritParams modulationSpectrum
#' @param amRes controls the width of window over which we look for local maxima
#' @return Returns a list with the frequency (Hz) and depth of amplitude
#' modulation (dB relative to global max, normally at 0 Hz).
#' @keywords internal
getAM = function(m,
amRange = c(10, 100),
amRes = NULL) {
if (is.null(amRes)) amRes = 0
colNames = abs(as.numeric(colnames(m)))
out = list(amMsFreq = NA, amMsPurity = NA)
# image(t(log(m)))
# fold around 0 and average AM across all FM bins
am = data.frame(freq = abs(colNames), amp = colSums(m))
am = am[order(am$freq), ]
# plot(am, type = 'l')
# average folded AM function cross pos/neg AM frequencies
# (b/c each point is doubled)
am_sm = am
i = 1
while(i < nrow(am_sm)) {
if (abs(am_sm$freq[i] - am_sm$freq[i + 1]) < .1) {
am_sm$amp[i] = (am_sm$amp[i] + am_sm$amp[i + 1]) / 2
am_sm$amp[i + 1] = NA
i = i + 2
} else {
i = i + 1
}
}
am_sm = na.omit(am_sm)
# plot(am, type = 'l')
# lines(am_sm, type = 'l', col = 'blue')
# find local maxima within amRange (note that we don't just throw away the
# rest of ms to improve the precision of smoothed AM function)
am_smRan = am_sm[am_sm$freq >= amRange[1] & am_sm$freq <= amRange[2], ]
# plot(am_smRan, type = 'l')
wl = max(3, round(amRes / (colNames[2] - colNames[1])))
# eg if amRes = 10, we look for a local maximum within ±5 Hz
temp = zoo::rollapply(
zoo::as.zoo(am_smRan$amp),
width = wl, # width 10 Hz, ie ±5 Hz
align = 'center',
function(x) {
middle = ceiling(length(x) / 2)
return(which.max(x) == middle)
})
idx = zoo::index(temp)[zoo::coredata(temp)]
if (length(idx) > 0) {
peaks = am_smRan[idx, ]
peaks = peaks[which.max(peaks$amp), ]
if (nrow(peaks) > 0) {
out$amMsFreq = peaks$freq
out$amMsPurity = log10(peaks$amp / max(am_sm$amp)) * 20
}
}
return(out)
}
#' Log-warp a modulation spectrum
#'
#' Internal soundgen function
#'
#' Log-warps a modulation spectrum along time dimension
#'
#' @param x a modulation spectrum: rows = FM, cols = AM
#' @keywords internal
#' @examples
#' a = matrix(1:44, ncol = 11)
#' colnames(a) = -5:5
#' soundgen:::logWarpMS(a, logWarp = 2)
logWarpMS = function(x, logWarp) {
X = as.numeric(colnames(x))
neg_col = which(X < 0)
zero_col = which(X == 0)
pos_col = which(X > 0)
m_left = logMatrix(x[, rev(neg_col)], base = logWarp)
# NB: flip the left half!
m_right = logMatrix(x[, pos_col], base = logWarp)
x_transf = cbind(m_left[, ncol(m_left):1], x[, zero_col, drop = FALSE], m_right)
return(x_transf)
}
Try the soundgen package in your browser
Any scripts or data that you put into this service are public.
soundgen documentation built on Sept. 29, 2023, 5:09 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions logWarpMS getAM averageMatrices getRough modulationSpectrumFragment plotMS .modulationSpectrum modulationSpectrum
Documented in averageMatrices getAM getRough logWarpMS modulationSpectrum .modulationSpectrum modulationSpectrumFragment plotMS
### MODULATION SPECTRUM
#' Modulation spectrum
#'
#' Produces a modulation spectrum of waveform(s) or audio file(s), with temporal
#' modulation along the X axis (Hz) and spectral modulation (1/KHz) along the Y
#' axis. A good visual analogy is decomposing the spectrogram into a sum of
#' ripples of various frequencies and directions. Roughness is calculated as the
#' proportion of energy / amplitude of the modulation spectrum within
#' \code{roughRange} of temporal modulation frequencies. The frequency of
#' amplitude modulation (amMsFreq, Hz) is calculated as the highest peak in the
#' smoothed AM function, and its purity (amMsPurity, dB) as the ratio of this
#' peak to the median AM over \code{amRange}. For relatively short and steady
#' sounds, set \code{amRes = NULL} and analyze the entire sound. For longer
#' sounds and when roughness or AM vary over time, set \code{amRes} to get
#' multiple measurements over time (see examples).
#'
#' Algorithm: prepare a spectrogram, take its logarithm (if \code{logSpec =
#' TRUE}), center, perform a 2D Fourier transform (see also
#' spectral::spec.fft()), take the upper half of the resulting symmetric matrix,
#' and raise it to \code{power}. The result is returned as \code{$original}. For
#' plotting purposes, the modulation matrix can be smoothed with Gaussian blur
#' (see \code{\link{gaussianSmooth2D}}) and log-warped (if \code{logWarp} is a
#' positive number). This processed modulation spectrum is returned as
#' \code{$processed}. If the audio is long enough, multiple windows are
#' analyzed, resulting in a vector of roughness values. For multiple inputs,
#' such as a list of waveforms or path to a folder with audio files, the
#' ensemble of modulation spectra can be interpolated to the same spectral and
#' temporal resolution and averaged (if \code{averageMS}).
#'
#' @seealso \code{\link{plotMS}} \code{\link{spectrogram}} \code{\link{analyze}}
#'
#' @references \itemize{
#' \item Singh, N. C., & Theunissen, F. E. (2003). Modulation spectra of
#' natural sounds and ethological theories of auditory processing. The Journal
#' of the Acoustical Society of America, 114(6), 3394-3411.
#' }
#' @inheritParams spectrogram
#' @inheritParams analyze
#' @param roughRange the range of temporal modulation frequencies that
#' constitute the "roughness" zone, Hz
#' @param amRange the range of temporal modulation frequencies that we are
#' interested in as "amplitude modulation" (AM), Hz
#' @param amRes target resolution of amplitude modulation, Hz. If \code{NULL},
#' the entire sound is analyzed at once, resulting in a single roughness value
#' (unless it is longer than \code{maxDur}, in which case it is analyzed in
#' chunks \code{maxDur} s long). If \code{amRes} is set, roughness is
#' calculated for windows \code{~1000/amRes} ms long (but at least 3 STFT
#' frames). \code{amRes} also affects the amount of smoothing when calculating
#' \code{amMsFreq} and \code{amMsPurity}
#' @param maxDur sounds longer than \code{maxDur} s are split into fragments,
#' and the modulation spectra of all fragments are averaged
#' @param logSpec if TRUE, the spectrogram is log-transformed prior to taking 2D
#' FFT
#' @param power raise modulation spectrum to this power (eg power = 2 for ^2, or
#' "power spectrum")
#' @param normalize if TRUE, the modulation spectrum of each analyzed fragment
#' \code{maxDur} in duration is separately normalized to have max = 1
#' @param returnMS if FALSE, only roughness is returned (much faster)
#' @param returnComplex if TRUE, returns a complex modulation spectrum (without
#' normalization and warping)
#' @param averageMS if TRUE, the modulation spectra of all inputs are averaged
#' into a single output; if FALSE, a separate MS is returned for each input
#' @param plot if TRUE, plots the modulation spectrum of each sound
#' @param savePlots if a valid path is specified, a plot is saved in this folder
#' (defaults to NA)
#' @param logWarpX,logWarpY numeric vector of length 2: c(sigma, base) of
#' pseudolog-warping the modulation spectrum, as in function
#' pseudo_log_trans() from the "scales" package
#' @param quantiles labeled contour values, \% (e.g., "50" marks regions that
#' contain 50\% of the sum total of the entire modulation spectrum)
#' @param kernelSize the size of Gaussian kernel used for smoothing (1 = no
#' smoothing)
#' @param kernelSD the SD of Gaussian kernel used for smoothing, relative to its
#' size
#' @param xlab,ylab,main,xlim,ylim graphical parameters
#' @param width,height,units,res parameters passed to
#' \code{\link[grDevices]{png}} if the plot is saved
#' @param ... other graphical parameters passed on to \code{filled.contour.mod}
#' and \code{\link[graphics]{contour}} (see \code{\link{spectrogram}})
#' @return Returns a list with the following components:
#' \itemize{
#' \item \code{$original} modulation spectrum prior to blurring and log-warping,
#' but after squaring if \code{power = TRUE}, a matrix of nonnegative values.
#' Rownames are spectral modulation frequencies (cycles/KHz), and colnames are
#' temporal modulation frequencies (Hz).
#' \item \code{$processed} modulation spectrum after blurring and log-warping
#' \item \code{$complex} untransformed complex modulation spectrum (returned
#' only if returnComplex = TRUE)
#' \item \code{$roughness} proportion of energy / amplitude of the modulation
#' spectrum within \code{roughRange} of temporal modulation frequencies, \% - a
#' vector if amRes is numeric and the sound is long enough, a single number
#' otherwise
#' \item \code{$amMsFreq} frequency of the highest peak, within \code{amRange}, of
#' the folded AM function (average AM across all FM bins for both negative and
#' positive AM frequencies), where a peak is a local maximum over \code{amRes}
#' Hz. Like \code{roughness}, \code{amMsFreq} and \code{amMsPurity} can be single
#' numbers or vectors, depending on whether the sound is analyzed as a whole or
#' in chunks
#' \item \code{$amMsPurity} ratio of the peak at amMsFreq to the median AM over
#' \code{amRange}, dB
#' \item \code{$summary} dataframe with summaries of roughness, amMsFreq, and
#' amMsPurity
#' }
#' @export
#' @examples
#' # White noise
#' ms = modulationSpectrum(runif(16000), samplingRate = 16000,
#' logSpec = FALSE, power = TRUE,
#' amRes = NULL) # analyze the entire sound, giving a single roughness value
#' str(ms)
#'
#' # Harmonic sound
#' s = soundgen(amFreq = 25, amDep = 50)
#' ms = modulationSpectrum(s, samplingRate = 16000, amRes = NULL)
#' ms[c('roughness', 'amMsFreq', 'amMsPurity')] # a single value for each
#' ms1 = modulationSpectrum(s, samplingRate = 16000, amRes = 5)
#' ms1[c('roughness', 'amMsFreq', 'amMsPurity')]
#' # measured over time (low values of amRes mean more precision, so we analyze
#' # longer segments and get fewer values per sound)
#'
#' # Embellish
#' ms = modulationSpectrum(s, samplingRate = 16000,
#' xlab = 'Temporal modulation, Hz', ylab = 'Spectral modulation, 1/KHz',
#' colorTheme = 'heat.colors', main = 'Modulation spectrum', lty = 3)
#'
#' \dontrun{
#' # A long sound with varying AM and a bit of chaos at the end
#' s_long = soundgen(sylLen = 3500, pitch = c(250, 320, 280),
#' amFreq = c(30, 55), amDep = c(20, 60, 40),
#' jitterDep = c(0, 0, 2))
#' playme(s_long)
#' ms = modulationSpectrum(s_long, 16000)
#' # plot AM over time
#' plot(x = seq(1, 1500, length.out = length(ms$amMsFreq)), y = ms$amMsFreq,
#' cex = 10^(ms$amMsPurity/20) * 10, xlab = 'Time, ms', ylab = 'AM frequency, Hz')
#' # plot roughness over time
#' spectrogram(s_long, 16000, ylim = c(0, 4),
#' extraContour = list(ms$roughness / max(ms$roughness) * 4000, col = 'blue'))
#'
#' # As with spectrograms, there is a tradeoff in time-frequency resolution
#' s = soundgen(pitch = 500, amFreq = 50, amDep = 100,
#' samplingRate = 44100, plot = TRUE)
#' # playme(s, samplingRate = 44100)
#' ms = modulationSpectrum(s, samplingRate = 44100,
#' windowLength = 50, step = 50, amRes = NULL) # poor temporal resolution
#' ms = modulationSpectrum(s, samplingRate = 44100,
#' windowLength = 5, step = 1, amRes = NULL) # poor frequency resolution
#' ms = modulationSpectrum(s, samplingRate = 44100,
#' windowLength = 15, step = 3, amRes = NULL) # a reasonable compromise
#'
#' # customize the plot
#' ms = modulationSpectrum(s, samplingRate = 44100,
#' windowLength = 15, overlap = 80, amRes = NULL,
#' kernelSize = 17, # more smoothing
#' xlim = c(-70, 70), ylim = c(0, 4), # zoom in on the central region
#' quantiles = c(.25, .5, .8), # customize contour lines
#' col = rev(rainbow(100)), # alternative palette
#' power = 2) # ^2
#' # Note the peaks at FM = 2/KHz (from "pitch = 500") and AM = 50 Hz (from
#' # "amFreq = 50")
#'
#' # Input can be a wav/mp3 file
#' ms = modulationSpectrum('~/Downloads/temp/16002_Faking_It_Large_clear.wav')
#'
#' # Input can be path to folder with audio files. Each file is processed
#' # separately, and the output can contain an MS per file...
#' ms1 = modulationSpectrum('~/Downloads/temp', kernelSize = 11,
#' plot = FALSE, averageMS = FALSE)
#' ms1$summary
#' names(ms1$original) # a separate MS per file
#' # ...or a single MS can be calculated:
#' ms2 = modulationSpectrum('~/Downloads/temp', kernelSize = 11,
#' plot = FALSE, averageMS = TRUE)
#' plotMS(ms2$original)
#' ms2$summary
#'
#' # Input can also be a list of waveforms (numeric vectors)
#' ss = vector('list', 10)
#' for (i in 1:length(ss)) {
#' ss[[i]] = soundgen(sylLen = runif(1, 100, 1000), temperature = .4,
#' pitch = runif(3, 400, 600))
#' }
#' # lapply(ss, playme)
#' # MS of the first sound
#' ms1 = modulationSpectrum(ss[[1]], samplingRate = 16000, scale = 1)
#' # average MS of all 10 sounds
#' ms2 = modulationSpectrum(ss, samplingRate = 16000, scale = 1, averageMS = TRUE, plot = FALSE)
#' plotMS(ms2$original)
#'
#' # A sound with ~3 syllables per second and only downsweeps in F0 contour
#' s = soundgen(nSyl = 8, sylLen = 200, pauseLen = 100, pitch = c(300, 200))
#' # playme(s)
#' ms = modulationSpectrum(s, samplingRate = 16000, maxDur = .5,
#' xlim = c(-25, 25), colorTheme = 'seewave',
#' power = 2)
#' # note the asymmetry b/c of downsweeps
#'
#' # "power = 2" returns squared modulation spectrum - note that this affects
#' # the roughness measure!
#' ms$roughness
#' # compare:
#' modulationSpectrum(s, samplingRate = 16000, maxDur = .5,
#' xlim = c(-25, 25), colorTheme = 'seewave',
#' power = 1)$roughness # much higher roughness
#'
#' # Plotting with or without log-warping the modulation spectrum:
#' ms = modulationSpectrum(soundgen(), samplingRate = 16000, plot = TRUE)
#' ms = modulationSpectrum(soundgen(), samplingRate = 16000,
#' logWarpX = c(2, 2), plot = TRUE)
#'
#' # logWarp and kernelSize have no effect on roughness
#' # because it is calculated before these transforms:
#' modulationSpectrum(s, samplingRate = 16000, logWarpX = c(1, 10))$roughness
#' modulationSpectrum(s, samplingRate = 16000, logWarpX = NA)$roughness
#' modulationSpectrum(s, samplingRate = 16000, kernelSize = 17)$roughness
#'
#' # Log-transform the spectrogram prior to 2D FFT (affects roughness):
#' modulationSpectrum(s, samplingRate = 16000, logSpec = FALSE)$roughness
#' modulationSpectrum(s, samplingRate = 16000, logSpec = TRUE)$roughness
#'
#' # Complex modulation spectrum with phase preserved
#' ms = modulationSpectrum(soundgen(), samplingRate = 16000,
#' returnComplex = TRUE)
#' plotMS(abs(ms$complex)) # note the symmetry
#' # compare:
#' plotMS(ms$original)
#' }
modulationSpectrum = function(
x,
samplingRate = NULL,
scale = NULL,
from = NULL,
to = NULL,
amRes = 5,
maxDur = 5,
logSpec = FALSE,
windowLength = 15,
step = NULL,
overlap = 80,
wn = 'hanning',
zp = 0,
power = 1,
normalize = TRUE,
roughRange = c(30, 150),
amRange = c(10, 200),
returnMS = TRUE,
returnComplex = FALSE,
summaryFun = c('mean', 'median', 'sd'),
averageMS = FALSE,
reportEvery = NULL,
cores = 1,
plot = TRUE,
savePlots = NULL,
logWarpX = NULL,
logWarpY = NULL,
quantiles = c(.5, .8, .9),
kernelSize = 5,
kernelSD = .5,
colorTheme = c('bw', 'seewave', 'heat.colors', '...')[1],
col = NULL,
main = NULL,
xlab = 'Hz',
ylab = '1/KHz',
xlim = NULL,
ylim = NULL,
width = 900,
height = 500,
units = 'px',
res = NA,
...
) {
## Prepare a list of arguments to pass to .modulationSpectrum()
myPars = c(as.list(environment()), list(...))
# exclude unnecessary args
myPars = myPars[!names(myPars) %in% c(
'x', 'samplingRate', 'scale', 'from', 'to', 'savePlots',
'reportEvery', 'cores', 'summaryFun', 'averageMS')]
# analyze
pa = processAudio(
x,
samplingRate = samplingRate,
scale = scale,
from = from,
to = to,
funToCall = '.modulationSpectrum',
myPars = myPars,
reportEvery = reportEvery,
cores = cores,
savePlots = savePlots
)
# htmlPlots
if (!is.null(pa$input$savePlots) && pa$input$n > 1) {
try(htmlPlots(pa$input, savePlots = savePlots, changesAudio = FALSE,
suffix = "MS", width = paste0(width, units)))
}
# prepare output
if (!is.null(summaryFun) && any(!is.na(summaryFun))) {
temp = vector('list', pa$input$n)
for (i in 1:pa$input$n) {
if (!pa$input$failed[i]) {
temp[[i]] = summarizeAnalyze(
data.frame(roughness = pa$result[[i]]$roughness,
amMsFreq = pa$result[[i]]$amMsFreq,
amMsPurity = pa$result[[i]]$amMsPurity),
summaryFun = summaryFun,
var_noSummary = NULL)
}
}
idx_failed = which(pa$input$failed)
idx_ok = which(!pa$input$failed)
if (length(idx_failed) > 0) {
idx_ok = which(!pa$input$failed)
if (length(idx_ok) > 0) {
filler = temp[[idx_ok[1]]] [1, ]
filler[1, ] = NA
} else {
stop('Failed to analyze any input')
}
for (i in idx_failed) temp[[i]] = filler
}
mysum_all = cbind(data.frame(file = pa$input$filenames_base),
do.call('rbind', temp))
} else {
mysum_all = NULL
}
for (i in idx_failed) pa$result[[i]] = list(
original = NULL, processed = NULL, complex = NULL,
roughness = NULL, amMsFreq = NULL, amMsPurity = NULL
)
original = processed = complex = roughness = amMsFreq = amMsPurity = NULL
# (otherwise note about no visible binding)
out_prep = c('original', 'processed', 'complex', 'roughness', 'amMsFreq', 'amMsPurity')
if (pa$input$n == 1) {
# unlist
for (op in out_prep) assign(noquote(op), pa$result[[1]] [[op]])
} else {
for (op in out_prep) assign(noquote(op), lapply(pa$result, function(x) x[[op]]))
}
# average MS across sounds
op_prep2 = c('original', 'processed') # roughness is not a matrix
if (returnComplex) {
op_prep2 = c(op_prep2, 'complex')
} else {
complex = NULL
}
if (averageMS & pa$input$n > 1) {
for (op in op_prep2) {
assign(
noquote(op),
averageMatrices(
get(op) [idx_ok],
rFun = 'max', # normally same samplingRate, but if not, upsample frequency resolution
cFun = 'median', # typical ncol (depends on sound dur)
reduceFun = '+') # internally divides by length, so actually becomes mean
)
}
}
invisible(list(original = original,
processed = processed,
complex = complex,
roughness = roughness,
amMsFreq = amMsFreq,
amMsPurity = amMsPurity,
summary = mysum_all))
}
#' Modulation spectrum per sound
#'
#' Internal soundgen function.
#' @inheritParams modulationSpectrum
#' @param audio a list returned by \code{readAudio}
#' @keywords internal
.modulationSpectrum = function(
audio,
amRes = 5,
maxDur = 5,
logSpec = FALSE,
windowLength = 15,
step = NULL,
overlap = 80,
wn = 'hanning',
zp = 0,
power = 1,
normalize = TRUE,
roughRange = c(30, 150),
amRange = c(10, 200),
returnMS = TRUE,
returnComplex = FALSE,
plot = TRUE,
savePlots = NULL,
logWarpX = NULL,
logWarpY = NULL,
quantiles = c(.5, .8, .9),
kernelSize = 5,
kernelSD = .5,
colorTheme = c('bw', 'seewave', 'heat.colors', '...')[1],
col = NULL,
main = NULL,
xlab = 'Hz',
ylab = '1/KHz',
xlim = NULL,
ylim = NULL,
width = 900,
height = 500,
units = 'px',
res = NA,
...
) {
# Re-set windowLength, step, and overlap so as to ensure that
# windowLength_points and step_points are not fractions
if (is.null(step)) step = windowLength * (1 - overlap / 100)
step_points = round(step / 1000 * audio$samplingRate)
step = step_points / audio$samplingRate * 1000
windowLength_points = round(windowLength / 1000 * audio$samplingRate)
windowLength = windowLength_points / audio$samplingRate * 1000
overlap = 100 * (1 - step_points / windowLength_points)
max_am = 1000 / step / 2
if (max_am < roughRange[1]) {
warning(paste(
'roughRange outside the analyzed range of temporal modulation frequencies;',
'increase overlap / decrease step to improve temporal resolution,',
'or else look for roughness in a lower range'))
}
if (is.numeric(amRes)) {
nFrames = max(3, ceiling(max_am / amRes * 2)) # min 3 spectrograms frames
} else {
nFrames = NULL
}
if (is.numeric(nFrames)) {
# split the input sound into chunks nFrames long
chunk_ms = windowLength + step * (nFrames - 1)
splitInto = max(1, ceiling(audio$duration * 1000 / chunk_ms))
if (chunk_ms > (audio$duration * 1000)) {
message(paste('The sound is too short to be analyzed with amRes =', amRes,
'Hz. Roughness is probably not measured correctly'))
chunk_ms = audio$duration
}
} else {
# split only those sounds that exceed maxDur
splitInto = max(1, ceiling(audio$duration / maxDur))
# so, for ex., if 2.1 times longer than maxDur, split into three
}
if (splitInto > 1) {
myseq = floor(seq(1, length(audio$sound), length.out = splitInto + 1))
myInput = vector('list', splitInto)
for (i in 1:splitInto) {
idx = myseq[i]:(myseq[i + 1])
myInput[[i]] = audio$sound[idx]
}
} else {
myInput = list(audio$sound)
}
# extract modulation spectrum per fragment
out = vector('list', splitInto)
if (returnComplex) {
out_complex = out
} else {
out_aggreg_complex = NULL
}
for (i in 1:splitInto) {
ms_i = modulationSpectrumFragment(
myInput[[i]],
samplingRate = audio$samplingRate,
windowLength = windowLength,
windowLength_points = windowLength_points,
step = step,
step_points = step_points,
wn = wn,
zp = zp,
logSpec = logSpec,
power = power,
normalize = normalize)
out[[i]] = ms_i$ms_half
if (returnComplex) {
out_complex[[i]] = ms_i$ms_complex
}
}
keep = which(unlist(lapply(out, function(x) !is.null(x))))
out = out[keep]
if (length(out) < 1) {
warning('The sound is too short or windowLength too long. Need at least 3 frames')
return(list('original' = NA,
'processed' = NA,
'complex' = NA,
'roughness' = NA,
'amMsFreq' = NA,
'amMsPurity' = NA))
}
# extract a measure of roughness
roughness = unlist(lapply(out, function(x)
getRough(x, roughRange, colNames = as.numeric(colnames(out[[1]])))))
# detect systematic amplitude modulation at the same frequency
am_list = lapply(out, function(x)
getAM(x, amRange = amRange, amRes = amRes))
amMsFreq = unlist(lapply(am_list, function(x) x$amMsFreq))
amMsPurity = unlist(lapply(am_list, function(x) x$amMsPurity))
# average modulation spectra across all sounds
if (!returnMS) {
result = list('original' = NULL,
'processed' = NULL,
'complex' = NULL,
'roughness' = roughness,
'amMsFreq' = amMsFreq,
'amMsPurity' = amMsPurity)
} else {
out_aggreg = averageMatrices(
out,
rFun = 'max', # normally same samplingRate, but if not, upsample frequency resolution
cFun = 'median', # typical ncol (depends on sound dur)
reduceFun = '+') # internally divides by length, so actually becomes mean
X = as.numeric(colnames(out_aggreg))
Y = as.numeric(rownames(out_aggreg))
# image(X, Y, t(log(out_aggreg)))
# prepare a separate summary of the complex ms
if (returnComplex) {
out_aggreg_complex = averageMatrices(
out_complex,
rFun = 'max',
cFun = 'median',
reduceFun = '+')
# image(t(log(abs(out_aggreg_complex))))
}
# smoothing / blurring
out_transf = gaussianSmooth2D(out_aggreg,
kernelSize = kernelSize,
kernelSD = kernelSD)
result = list('original' = out_aggreg,
'processed' = out_transf,
'complex' = out_aggreg_complex,
'roughness' = roughness,
'amMsFreq' = amMsFreq,
'amMsPurity' = amMsPurity)
} # end of if (returnMS)
# PLOTTING
if (is.character(audio$savePlots)) {
plot = TRUE
png(filename = paste0(audio$savePlots, audio$filename_noExt, "_MS.png"),
width = width, height = height, units = units, res = res)
}
if (plot) {
plotMS(ms = out_transf,
X = X,
Y = Y,
audio = audio,
colorTheme = colorTheme,
col = col,
logWarpX = logWarpX,
logWarpY = logWarpY,
xlab = xlab, ylab = ylab,
main = main,
xlim = xlim, ylim = ylim,
quantiles = quantiles,
...
)
if (is.character(audio$savePlots)) dev.off()
}
invisible(result)
}
#' Plot modulation spectrum
#'
#' Plots a single modulation spectrum returned by
#' \code{\link{modulationSpectrum}}. The result is the same as the plot produced
#' by \code{\link{modulationSpectrum}}, but calling \code{plotMS} is handy for
#' processed modulation spectra - for instance, for plotting the difference
#' between the modulation spectra of two sounds or groups of sounds.
#'
#' @inheritParams modulationSpectrum
#' @param ms modulation spectrum - a matrix with temporal modulation in columns
#' and spectral modulation in rows, as returned by
#' \code{\link{modulationSpectrum}}
#' @param X,Y rownames and colnames of \code{ms}, respectively
#' @param audio (internal) a list of audio attributes
#' @export
#' @examples
#' ms1 = modulationSpectrum(runif(4000), samplingRate = 16000, plot = TRUE)
#' plotMS(ms1$processed) # identical to above
#'
#' # compare two modulation spectra
#' ms2 = modulationSpectrum(soundgen(sylLen = 100, addSilence = 0),
#' samplingRate = 16000)
#' # ensure the two matrices have the same dimensions
#' ms2_resized = soundgen:::interpolMatrix(ms2$original,
#' nr = nrow(ms1$original), nc = ncol(ms1$original))
#' # plot the difference
#' plotMS(log(ms1$original / ms2_resized), quantile = NULL,
#' col = colorRampPalette(c('blue', 'yellow')) (50))
plotMS = function(
ms,
X = NULL,
Y = NULL,
quantiles = c(.5, .8, .9),
colorTheme = c('bw', 'seewave', 'heat.colors', '...')[1],
col = NULL,
logWarpX = NULL,
logWarpY = NULL,
main = NULL,
xlab = 'Hz',
ylab = '1/KHz',
xlim = NULL,
ylim = NULL,
audio = NULL,
...
) {
if (!is.null(col)) colorTheme = NULL
if (!is.null(colorTheme)) {
color.palette = switchColorTheme(colorTheme)
} else {
color.palette = NULL
}
if (is.null(X)) X = as.numeric(colnames(ms))
if (is.null(Y)) Y = as.numeric(rownames(ms))
if (is.null(xlim)) xlim = c(X[1], -X[1])
if (is.null(ylim)) ylim = range(Y)
if (is.null(main) & !is.null(audio$filename_noExt)) {
if (audio$filename_noExt == 'sound') {
main = ''
} else {
main = audio$filename_noExt
}
}
# plot with filled.contour.mod
ms = zeroOne(ms)
X1 = X
Y1 = Y
if (is.numeric(logWarpX) | is.numeric(logWarpY)) {
if (is.numeric(logWarpX)) {
X1 = pseudoLog(X, sigma = logWarpX[1], base = logWarpX[2])
# lab_x = pretty(X, n = 9, min.n = 7)
# at_x = pseudoLog(lab_x, sigma = logWarpX[1], base = logWarpX[2])
at_x = pretty(X1, n = 7)
lab_x = round(pseudoLog_undo(at_x, sigma = logWarpX[1], base = logWarpX[2]))
} else {
lab_x = at_x = pretty(X)
}
if (is.numeric(logWarpY)) {
Y1 = pseudoLog(Y, sigma = logWarpY[1], base = logWarpY[2])
lab_y = pretty(Y)
at_y = pseudoLog(lab_y, sigma = logWarpY[2], base = logWarpY[2])
} else {
lab_y = at_y = pretty(Y)
}
filled.contour.mod(X1, Y1, t(ms),
levels = seq(0, 1, length = 100),
color.palette = color.palette,
col = col,
xlab = xlab, ylab = ylab,
bty = 'n',
main = main, xaxt = 'n', yaxt = 'n',
# xlim = xlim, ylim = ylim,
...)
# add nicely labeled axes
axis(1, at = at_x, labels = lab_x, ...) # xpd = TRUE to extend beyond the plot
axis(2, at = at_y, labels = lab_y, ...)
lbl_Hz_pos = pretty(Y1)
lbls_Hz = round(1000 / lbl_Hz_pos)
axis(4, at = lbl_Hz_pos, labels = lbls_Hz, ...)
} else {
filled.contour.mod(X1, Y1, t(ms),
levels = seq(0, 1, length = 100),
color.palette = color.palette,
col = col,
xlab = xlab, ylab = ylab,
bty = 'n',
main = main,
# xlim = xlim, ylim = ylim,
...)
lbls = round(1000 / pretty(Y1))
lbl_pos = 1000 / lbls
axis(4, at = lbl_pos, labels = lbls, ...)
}
abline(v = 0, lty = 3)
# add contours
if (is.numeric(quantiles)) {
# qntls = quantile(out_aggreg, probs = quantiles) # could try HDI instead
qntls = pDistr(as.numeric(ms), quantiles = quantiles)
par(new = TRUE)
contour(x = X1, y = Y1, z = t(ms),
levels = qntls, labels = quantiles * 100,
xaxs = 'i', yaxs = 'i',
axes = FALSE, frame.plot = FALSE,
# xlim = xlim, ylim = ylim,
...)
par(new = FALSE)
}
}
#' Modulation spectrum per fragment
#'
#' Internal soundgen function.
#' @inheritParams modulationSpectrum
#' @param sound numeric vector
#' @keywords internal
#' @examples
#' s = soundgen(amFreq = 25, amDep = 100)
#' ms = soundgen:::modulationSpectrumFragment(s, 16000,
#' windowLength = 50, windowLength_points = .05 * 16000,
#' step = 5, step_points = .005 * 16000)
#' image(as.numeric(colnames(ms$ms_half)), as.numeric(rownames(ms$ms_half)),
#' t(log(ms$ms_half)))
modulationSpectrumFragment = function(sound,
samplingRate,
windowLength,
windowLength_points,
step,
step_points,
wn = 'hanning',
zp = 0,
logSpec = FALSE,
power = 1,
normalize = TRUE) {
# Calling stdft is ~80 times faster than going through spectrogram (!)
step_seq = seq(1, length(sound) + 1 - windowLength_points, step_points)
# print(length(step_seq))
if (length(step_seq) < 3) return(NULL)
s1 = seewave::stdft(
wave = as.matrix(sound),
wn = wn,
wl = windowLength_points,
f = samplingRate,
zp = zp,
step = step_seq,
scale = FALSE,
norm = FALSE,
complex = FALSE
)
# image(t(log(s1))); s1[1:3, 1:3]; dim(s1); range(s1)
# log-transform amplitudes
if (logSpec) {
positives = which(s1 > 0)
nonpositives = which(s1 <= 0)
s1[positives] = log(s1[positives])
if (length(positives) > 0 & length(nonpositives) > 0) {
s1[nonpositives] = min(s1[positives])
}
s1 = s1 - min(s1) + 1e-16 # positive
# image(t(log(s1)))
}
# 2D fft
ms_complex = specToMS(s1, windowLength = windowLength, step = step)
ms = abs(ms_complex)
# image(as.numeric(colnames(ms)), as.numeric(rownames(ms)), t(log(ms)))
symAxis = floor(nrow(ms) / 2) + 1
# ms[(symAxis - 2) : (symAxis + 2), 1:2]
ms_half = ms[symAxis:nrow(ms),, drop = FALSE] # take only the upper half (always even)
# power
if (is.numeric(power) && power != 1) ms_half = ms_half ^ power
# normalize
if (normalize && any(s1 != 0)) {
ms_half = ms_half - min(ms_half)
ms_half = ms_half / max(ms_half)
}
# image(x = as.numeric(colnames(ms_half)), z = t(log(ms_half)))
return(list(
ms_half = ms_half,
ms_complex = ms_complex
))
}
#' Get roughness
#'
#' Internal soundgen function
#'
#' Helper function for calculating roughness - the proportion of energy /
#' amplitude in the roughness range
#' @param m numeric matrix of non-negative values with colnames giving temporal
#' modulation frequency
#' @inheritParams modulationSpectrum
#' @param colNames numeric vector with frequencies corresponding to columns
#' @return Returns roughness in percent.
#' @keywords internal
#' @examples
#' m = modulationSpectrum(soundgen(), samplingRate = 16000)$original
#' soundgen:::getRough(m, roughRange = c(30, Inf))
getRough = function(m, roughRange, colNames = NULL) {
if (is.null(colNames)) colNames = abs(as.numeric(colnames(m)))
rough_cols = which(colNames > roughRange[1] &
colNames < roughRange[2])
if (length(rough_cols) > 0) {
roughness = sum(m[, rough_cols]) / sum(m) * 100
} else {
roughness = 0
}
# # alternatively / in addition, could try gaussian filter when calculating roughness
# ampl_filter = dnorm(colNames, mean = roughMean, sd = roughSD)
# ampl_filter = ampl_filter / sum(ampl_filter) # to pdf
# # plot(ampl_filter, type = 'l')
# roughness = sum(colSums(m) * ampl_filter) / sum(m) * 100
return(roughness)
}
#' Average matrices
#'
#' Internal soundgen function.
#'
#' Takes a list of matrices (normally modulation spectra), interpolates them to
#' have the same size, and then reduces them (eg takes the average).
#' @param mat_list a list of matrices to aggregate (eg spectrograms or
#' modulation spectra)
#' @param rFun, cFun functions used to determine the number of rows and columns
#' in the result
#' @param reduceFun function used to aggregate
#' @keywords internal
#' @examples
#' mat_list = list(
#' matrix(1:30, nrow = 5),
#' matrix(80:17, nrow = 8)
#' )
#' soundgen:::averageMatrices(mat_list)
#' soundgen:::averageMatrices(mat_list, cFun = 'max', reduceFun = '*')
averageMatrices = function(mat_list,
rFun = 'max',
cFun = 'median',
reduceFun = '+') {
# normally same samplingRate, but in case not, upsample frequency resolution
nr = round(do.call(rFun, list(unlist(lapply(mat_list, nrow)))))
# take typical ncol (depends on sound dur)
nc = round(do.call(cFun, list(unlist(lapply(mat_list, ncol)))))
mat_list_sameDim = lapply(mat_list, function(x) interpolMatrix(x, nr = nr, nc = nc))
agg = Reduce(reduceFun, mat_list_sameDim) / length(mat_list)
return(agg)
}
#' Get amplitude modulation
#'
#' Internal soundgen function
#'
#' Helper function for calculating amplitude modulation based on the modulation
#' spectrum. Algorithm: averages AM across all FM bins in the positive half of
#' the modulation spectrum and looks for a peak in the specified AM frequency
#' range.
#' @param m numeric matrix of non-negative values with colnames giving temporal
#' modulation frequency
#' @inheritParams modulationSpectrum
#' @param amRes controls the width of window over which we look for local maxima
#' @return Returns a list with the frequency (Hz) and depth of amplitude
#' modulation (dB relative to global max, normally at 0 Hz).
#' @keywords internal
getAM = function(m,
amRange = c(10, 100),
amRes = NULL) {
if (is.null(amRes)) amRes = 0
colNames = abs(as.numeric(colnames(m)))
out = list(amMsFreq = NA, amMsPurity = NA)
# image(t(log(m)))
# fold around 0 and average AM across all FM bins
am = data.frame(freq = abs(colNames), amp = colSums(m))
am = am[order(am$freq), ]
# plot(am, type = 'l')
# average folded AM function cross pos/neg AM frequencies
# (b/c each point is doubled)
am_sm = am
i = 1
while(i < nrow(am_sm)) {
if (abs(am_sm$freq[i] - am_sm$freq[i + 1]) < .1) {
am_sm$amp[i] = (am_sm$amp[i] + am_sm$amp[i + 1]) / 2
am_sm$amp[i + 1] = NA
i = i + 2
} else {
i = i + 1
}
}
am_sm = na.omit(am_sm)
# plot(am, type = 'l')
# lines(am_sm, type = 'l', col = 'blue')
# find local maxima within amRange (note that we don't just throw away the
# rest of ms to improve the precision of smoothed AM function)
am_smRan = am_sm[am_sm$freq >= amRange[1] & am_sm$freq <= amRange[2], ]
# plot(am_smRan, type = 'l')
wl = max(3, round(amRes / (colNames[2] - colNames[1])))
# eg if amRes = 10, we look for a local maximum within ±5 Hz
temp = zoo::rollapply(
zoo::as.zoo(am_smRan$amp),
width = wl, # width 10 Hz, ie ±5 Hz
align = 'center',
function(x) {
middle = ceiling(length(x) / 2)
return(which.max(x) == middle)
})
idx = zoo::index(temp)[zoo::coredata(temp)]
if (length(idx) > 0) {
peaks = am_smRan[idx, ]
peaks = peaks[which.max(peaks$amp), ]
if (nrow(peaks) > 0) {
out$amMsFreq = peaks$freq
out$amMsPurity = log10(peaks$amp / max(am_sm$amp)) * 20
}
}
return(out)
}
#' Log-warp a modulation spectrum
#'
#' Internal soundgen function
#'
#' Log-warps a modulation spectrum along time dimension
#'
#' @param x a modulation spectrum: rows = FM, cols = AM
#' @keywords internal
#' @examples
#' a = matrix(1:44, ncol = 11)
#' colnames(a) = -5:5
#' soundgen:::logWarpMS(a, logWarp = 2)
logWarpMS = function(x, logWarp) {
X = as.numeric(colnames(x))
neg_col = which(X < 0)
zero_col = which(X == 0)
pos_col = which(X > 0)
m_left = logMatrix(x[, rev(neg_col)], base = logWarp)
# NB: flip the left half!
m_right = logMatrix(x[, pos_col], base = logWarp)
x_transf = cbind(m_left[, ncol(m_left):1], x[, zero_col, drop = FALSE], m_right)
return(x_transf)
}
Try the soundgen package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.