getSurprisal: Get surprisal

In soundgen: Sound Synthesis and Acoustic Analysis

Get surprisal

Description

Tracks the unpredictability of spectro-temporal changes in a sound over time, returning continuous contours of Shannon surprisal ($info), Bayesian surprise ($kl for Kullback-Leibler divergence), and autocorrelation-based surprisal ($surprisal). This is an attempt to track auditory salience over time - that is, to identify parts of a sound that are likely to involuntarily attract the listener's attention.

Usage

getSurprisal(
  x,
  samplingRate = NULL,
  scale = NULL,
  from = NULL,
  to = NULL,
  winSurp = 2000,
  specFun = "audSpec",
  specFun_pars = list(),
  logSpec = TRUE,
  method = c("acf", "none"),
  sameLagAllFreqs = FALSE,
  weightByAmpl = TRUE,
  weightByPrecision = TRUE,
  onlyPeakAutocor = TRUE,
  rescale = FALSE,
  minProb = 1e-12,
  summaryFun = "mean",
  output = c("surprisal", "info", "kl"),
  reportEvery = NULL,
  cores = 1,
  plot = TRUE,
  savePlots = FALSE,
  embed = FALSE,
  osc = c("linear", "dB", "none"),
  heights = c(3, 1),
  ylim = NULL,
  maxPoints = c(1e+05, 5e+05),
  colorTheme = "bw",
  col = NULL,
  extraContour = list(col = "blue", lwd = 3, lty = 1),
  xlab = NULL,
  ylab = NULL,
  xaxp = NULL,
  mar = c(5.1, 4.1, 4.1, 2),
  main = NULL,
  grid = NULL,
  width = 900,
  height = 500,
  units = "px",
  res = NA,
  ...
)

Arguments

x	path to a folder, one or more wav or mp3 files c('file1.wav', 'file2.mp3'), Wave object, numeric vector, or a list of Wave objects or numeric vectors
samplingRate	sampling rate of x (only needed if x is a numeric vector)
scale	maximum possible amplitude of input, used to normalize the input vector (only needed if x is a numeric vector)
from, to	if NULL (default), analyzes the whole sound, otherwise from...to (s)
winSurp	surprisal analysis window, ms. Inf means "from sound onset"; windows shorter than 3 frames produce NA
specFun	the function used to extract a spectrogram-like feature matrix. Can be a string or a custom function that takes audio (numeric vector) as the first argument and returns a spectrogram-like matrix with time in columns and features in rows (see examples). A precomputed spectrogram-like matrix is also accepted (features in rows, time in columns [ms], numeric rownames for plotting). Supported strings: stft_simple 'stft' / STFT / 'stft_simple' (amplitude spectrogram) Parameters in specFun_pars: samplingRate (optional if frequency and time labels are not needed), wl (samples), step (samples), wn, zp, padWithSilence. spectrogram 'spectrogram' (amplitude spectrogram - a wrapper around stft_simple with more options). Parameters in specFun_pars: see spectrogram. powspec 'powerspec' (power spectrogram with tuneR). Parameters in specFun_pars: wintime (s) or windowLength (ms), steptime (s) or step (ms), dither. melfcc 'melspec' (mel-spectrogram), 'mfcc' / 'melfcc' (MFCCs). Parameters in specFun_pars: windowLength (ms), step (ms), nbands, maxfreq (Hz), MFCC (integer vector). audSpectrogram 'audSpectrogram' / 'audSpec' (auditory spectrogram). Parameters in specFun_pars: see audSpectrogram. getRMS 'getRMS' / 'rms' (RMS amplitude envelope). Parameters in specFun_pars: see getRMS. getEnv 'getEnv' / 'env' (various smoothed envelopes: RMS, analytical, peak, etc.). Parameters in specFun_pars: see getEnv.
specFun_pars	a list of parameters passed to specFun. Defaults for audSpectrogram, list(yScale = 'ERB', nFilters_oct = 6, step = 15, minFreq = 60); for melspec, windowLength = 20, step = 20, maxfreq = NULL, nbands = 128, MFCC = 2:13; for spectrogram, windowLength = 20, step = 20; for env, windowLength = 40, step = 20 converted to samples; for rms, windowLength = 40, step = 20. For convenience when using melspec, windowLength and step in ms are converted to wintime and steptime in seconds if no explicit values in seconds are supplied.
logSpec	if TRUE, the output of specFun is log-transformed prior to calculating surprisal, offsetting as needed to avoid non-positive values
method	affects $surprisal and $bestLag_mat only; has no effect on $info and $kl. acf = change in autocorrelation at the previously best lag after adding the final point; none = do not calculate $surprisal; $surprisal and $bestLag_mat are then NA.
sameLagAllFreqs	only for method = 'acf'. If TRUE, the bestLag is calculated by averaging the ACFs of all channels, and the same bestLag is used to calculate the surprisal in each frequency channel (we expect the same "rhythm" for all frequencies). If FALSE, the bestLag is calculated separately for each frequency channel (we can track different "rhythms" at different frequencies).
weightByAmpl	f TRUE, ACF averaging (when sameLagAllFreqs = TRUE) and aggregation of all per-channel contours (surprisal, info, infoW, kl, klW) are weighted by the maximum non-negative value per frequency channel in the current analysis window. For log-transformed, processed, or custom feature matrices, these weights are feature maxima, not necessarily physical amplitude.
weightByPrecision	if TRUE, ACF-based surprisal is weighted by the current autocorrelation, so deviations from a previous pattern are more surprising if this pattern is strong.
onlyPeakAutocor	if TRUE, only peaks of ACFs are considered (so bestLag can never be 1, and the first change after a string of static values results in surprisal = NA).
rescale	if TRUE, aggregated surprisal is normalized from (-Inf, Inf) to approximately (-1, 1) using tanh(surprisal / 2). surprisal_mat is unaffected.
minProb	minimum probability used to cap Shannon surprisal: info and infoW cannot exceed -log(minProb). Also used as a lower bound for the KL divergence before logging in kl and klW.
summaryFun	functions used to summarize each acoustic characteristic, eg c('mean', 'sd'); user-defined functions are fine (see examples); NAs are omitted automatically for mean/median/sd/min/max/range/sum, otherwise take care of NAs yourself
output	what to return, options: 'surprisal', 'loudness', 'dLoudness', 'surprisalLoudness', 'surprisal_mat', 'bestLag_mat', 'info', 'info_mat', 'infoW', 'infoW_mat', 'kl', 'kl_mat', 'klW', 'klW_mat', 'spectrogram', 'all' (see the Return section)
reportEvery	when processing multiple inputs, report estimated time left every reportEvery iterations (NULL = default, NA = don't report); see reportTime
cores	number of cores for parallel processing
plot	If TRUE, plots the feature matrix and the surprisal contour. For specFun = 'env', the surprisal contour is overlaid on an oscillogram.
savePlots	if TRUE, creates a subdirectory in the input directory (if input is a file or folder) or in the working directory (if input is a vector etc), named after the function (eg "spectrogram/"). All plots and audio files (if any) are saved in this new directory. If there are multiple inputs, an html notebook is also created for easy viewing and listening
embed	if TRUE and savePlots is set and there are multiple inputs, all saved images and audio (if any) are embedded in the exported html notebook for easy sharing; if FALSE, the html file links to separate images and audio files (but separate files are still saved). NB: for this to work, package "base64enc" must be installed
osc	"none" = no oscillogram; "linear" = on the original scale; "dB" = in decibels
heights	a vector of length two specifying the relative height of the spectrogram and the oscillogram (including time axes labels)
ylim	frequency range to plot, kHz (defaults to 0 to Nyquist frequency). NB: still in kHz, even if yScale = bark, mel, or ERB
maxPoints	the maximum number of "pixels" in the oscillogram (if any) and spectrogram; good for quickly plotting long audio files; defaults to c(1e5, 5e5); does not affect reassigned spectrograms
colorTheme	black and white ('bw'), as in seewave package ('seewave'), matlab-type palette ('matlab'), or any palette from palette such as 'heat.colors', 'cm.colors', etc
col	actual colors, eg rev(rainbow(100)) - see ?hcl.colors for colors in base R (overrides colorTheme)
extraContour	a vector of arbitrary length scaled in Hz (regardless of yScale, but nonlinear yScale also warps the contour) that will be plotted over the spectrogram (eg pitch contour); can also be a list with extra graphical parameters such as lwd, col, warp (FALSE = plot as is, TRUE = warp to conform to nonlinear yScale), etc. (see examples)
xlab, ylab, main, mar, xaxp	graphical parameters for plotting
grid	if numeric, adds n = grid dotted lines per kHz and the same number of lines along the time axis
width, height, units, res	graphical parameters for saving plots passed to png
...	other graphical parameters

Details

Algorithm: the sound is transformed into some spectrogram-like representation (e.g., an auditory spectrogram, a mel-warped STFT spectrogram, etc.) or an RMS amplitude envelope. Using just the envelope is very fast, but then we discard all spectral information. For each frequency channel, a sliding window is analyzed to compare the actually observed final value with its expected value. The resulting per-channel surprisal contours are aggregated by taking their mean - optionally, weighted by the maximum amplitude of each frequency channel across the analysis window. Because increases in loudness are known to be important predictors of auditory salience, loudness per frame is also returned, as well as the product of its positive changes and surprisal.

Value

A list with two top-level elements: $detailed and $summary.

$detailed contains per-frame statistics selected with the output argument. If multiple sounds are analyzed, $detailed is a list of per-sound lists.

$summary contains per-file summaries of a fixed set of contours: loudness, surprisal, surprisalLoudness, info, infoW, kl, and klW. These are summarized even if not all of them are included in output. If summaryFun is NULL, $summary is NULL.

Available measures:

surprisal

Aggregated surprisal contour: change in autocorrelation. Values are averaged across frequency channels, optionally weighted by channel amplitude. Positive values mean "an unexpected change", and negative values mean "a change that confirms the expectations, making signal periodicity more certain." If rescale = TRUE, the aggregated contour is transformed with tanh(surprisal / 2) to approximately (-1, 1).

loudness

Subjective loudness in sone, as per getLoudness, resampled to match the number of surprisal frames.

dLoudness

First temporal derivative of max-normalized loudness: diff(c(0, loudness / max(loudness))).

surprisalLoudness

Product of the positive parts of surprisal and dLoudness: max(surprisal, 0) * max(dLoudness, 0). This contour emphasizes surprising events that coincide with increases in loudness. If rescale = TRUE, this uses the rescaled surprisal contour.

surprisal_mat

Matrix of per-channel surprisal values before aggregation across frequency channels (frequency channels in rows, time frames in columns). Column names are time in ms.

bestLag_mat

Matrix of the autocorrelation lag used to calculate ACF-surprisal in each time-frequency bin, in seconds. NA where no lag was available or applicable, static analysis windows, or when onlyPeakAutocor = TRUE and no ACF peak was found. If sameLagAllFreqs = TRUE, the same lag is used for all non-static frequency channels; static channels still return NA.

info

Shannon surprisal contour, calculated as -log(rho), where rho is the Gaussian density at the next observation normalized by the maximum Gaussian density. Values are capped at -log(minProb).

info_mat

Matrix of per-channel Shannon surprisal values corresponding to info.

infoW

Windowed Shannon surprisal: same as info, but using weighted means and standard deviations from a half-Gaussian taper that prioritizes more recent observations.

infoW_mat

Matrix of per-channel windowed Shannon surprisal values corresponding to infoW.

kl

Bayesian log-surprisal: natural logarithm of the Kullback-Leibler divergence between the Gaussian distributions before and after observing the next data point, with a window-length correction added as 2 * log(n). To avoid -Inf, the KL divergence is floored at minProb before taking the logarithm. Values can therefore be negative.

kl_mat

Matrix of per-channel Bayesian log-surprisal values corresponding to kl.

klW

Windowed Bayesian log-surprisal: same as kl, but using weighted means and variances from a half-Gaussian taper that prioritizes more recent observations. The full window length n is still used for the 2 * log(n) correction.

klW_mat

Matrix of per-channel windowed Bayesian log-surprisal values corresponding to klW.

spectrogram

The spectrogram-like feature matrix actually analyzed (frequency channels or features in rows, time frames in columns), after any requested preprocessing such as log-transformation. Column names are time in ms. If specFun = 'env', this is a one-row matrix containing the RMS envelope, possibly log-transformed if logSpec = TRUE.

References

• Anikin, A. (2026) Measuring surprisal in sound sequences. Behavior Research Methods. doi: 10.3758/s13428-026-03153-3.

Examples

# A quick example
data('speechEx', package = 'soundgen')
surp = getSurprisal(speechEx, from = 0.5, to = 1)
surp

## Not run: 
# A few more meaningful examples

## Example 1: a temporal deviant
s0 = soundgen(nSyl = 8, sylLen = 150,
              pauseLen = c(rep(200, 7), 450), pitch = c(200, 150),
              temperature = .05, plot = FALSE)
sound = c(rep(0, 4000),
          addVectors(rnorm(16000 * 3.5, 0, .02), s0, insertionPoint = 4000),
          rep(0, 200))
spectrogram(sound, 16000, yScale = 'ERB')

# long window (Inf = from the beginning)
surp = getSurprisal(sound, 16000, winSurp = Inf, output = 'all')
plot(sound, type = 'l')
surp_cont = surp$detailed$surprisal
lines(seq(0, length(sound), length.out = length(surp_cont)),
  surp_cont / max(surp_cont, na.rm = TRUE), col = 'blue', lwd = 2)

# Which frequency-time bins are surprising?
filled.contour(x = as.numeric(colnames(surp$detailed$surprisal_mat)) / 1000,
               y = as.numeric(rownames(surp$detailed$surprisal_mat)),
               z = t(surp$detailed$surprisal_mat),
               xlab = 'Time, s',
               ylab = 'Frequency, kHz')
# Best lag (periodicity) over time
hist(surp$detailed$bestLag_mat, xlab = 'Period, s')
abline(v = .35, lty = 3, lwd = 3, col = 'blue')  # true period = 350 ms

# just use the amplitude envelope instead of an auditory spectrogram
surp = getSurprisal(sound, 16000, winSurp = Inf, specFun = 'env')

# increase spectral and temporal resolution (can be slow)
surp = getSurprisal(sound, 16000, winSurp = 2000,
  specFun_pars = list(nFilters = 50, step = 10,
  yScale = 'bark', bandwidth = 1/4), output = 'all')

# weight by increase in loudness
spectrogram(sound, 16000, extraContour = surp$detailed$surprisalLoudness /
  max(surp$detailed$surprisalLoudness, na.rm = TRUE) * 8000)

par(mfrow = c(3, 1))
plot(surp$detailed$surprisal, type = 'l', xlab = '',
  ylab = '', main = 'surprisal')
abline(h = 0, lty = 2)
plot(surp$detailed$dLoudness, type = 'l', xlab = '',
  ylab = '', main = 'd-loudness')
abline(h = 0, lty = 2)
plot(surp$detailed$surprisalLoudness, type = 'l', xlab = '',
  ylab = '', main = 'surprisal * d-loudness')
par(mfrow = c(1, 1))

# short window = amnesia (every new sound is surprising)
getSurprisal(sound, 16000, winSurp = 300)

# add bells and whistles
surp = getSurprisal(sound, samplingRate = 16000,
  osc = 'dB',  # plot oscillogram in dB
  heights = c(2, 1),  # spectro/osc height ratio
  # colorTheme = 'heat.colors',  # pick color theme...
  col = rev(hcl.colors(30, palette = 'Viridis')),  # ...or specify the colors
  cex.lab = .75, cex.axis = .75,  # text size and other base graphics pars
  ylim = c(0, 5),  # always in kHz
  main = 'Audiogram with surprisal contour', # title
  extraContour = list(col = 'blue', lty = 2, lwd = 2)
  # + axis labels, etc
)

## Example 2: a spectral deviant
s1 = soundgen(
  nSyl = 11, sylLen = 150, invalidArgAction = 'ignore',
  formants = NULL, lipRad = 0,  # so all syls have the same envelope
  pauseLen = 90, pitch = c(1000, 750), rolloff = -20,
  pitchGlobal = c(rep(0, 5), 18, rep(0, 5)),
  temperature = .01, pitchCeiling = 7000,
  plot = TRUE, windowLength = 35)
surp = getSurprisal(s1, 16000, winSurp = 1500, output = 'all')
filled.contour(x = as.numeric(colnames(surp$detailed$surprisal_mat)) / 1000,
               y = as.numeric(rownames(surp$detailed$surprisal_mat)),
               z = t(surp$detailed$surprisal_mat),
               xlab = 'Time, s',
               ylab = 'Frequency, kHz')
# deviant surprising both at 1 kHz (expected tone omitted) and at the new freq
surp = getSurprisal(s1, 16000, winSurp = 1500,
  specFun = 'env')  # doesn't work - need spectral info

## Example 3: different rhythms in different frequency bins
s6_1 = soundgen(nSyl = 23, sylLen = 100, pauseLen = 50, pitch = 1200,
  rolloffExact = 1, invalidArgAction = 'ignore', plot = TRUE)
s6_2 = soundgen(nSyl = 10, sylLen = 250, pauseLen = 100, pitch = 400,
  rolloffExact = 1, invalidArgAction = 'ignore', plot = TRUE)
s6_3 = soundgen(nSyl = 5, sylLen = 400, pauseLen = 200, pitch = 3400,
  rolloffExact = 1, invalidArgAction = 'ignore', plot = TRUE)
s6 = addVectors(s6_1, s6_2)
s6 = addVectors(s6, s6_3)

surp = getSurprisal(s6, 16000, winSurp = Inf, sameLagAllFreqs = TRUE,
  specFun_pars = list(nFilters = 32), output = 'all')
surp = getSurprisal(s6, 16000, winSurp = Inf, sameLagAllFreqs = FALSE,
  specFun_pars = list(nFilters = 32), output = 'all')  # learns all 3 rhythms
filled.contour(x = as.numeric(colnames(surp$detailed$surprisal_mat)) / 1000,
               y = as.numeric(rownames(surp$detailed$surprisal_mat)),
               z = t(surp$detailed$surprisal_mat),
               xlab = 'Time, s',
               ylab = 'Frequency, kHz')

## Example 4: different time scales
s8 = soundgen(nSyl = 4, sylLen = 75, pauseLen = 50)
s8 = rep(c(s8, rep(0, 2000)), 8)
getSurprisal(s8, 16000, specFun = 'env', winSurp = Inf)
# ACF picks up first the fast rhythm, then after a few cycles switches to
# the slow rhythm

# Custom input: produce a nice spectrogram first, then use it as input
sp = spectrogram(s0, 16000, windowLength = 10, step = 10, contrast = .3,
  output = 'processed')  # return the modified spectrogram
colnames(sp) = as.numeric(colnames(sp)) / 1000  # convert ms to s
getSurprisal(s0, 16000, specFun = sp, logSpec = FALSE)

# Custom input: use acoustic features returned by analyze()
an = analyze(sound, 16000, windowLength = 20, novelty = NULL)
feature_mat = t(an$detailed[, 4:ncol(an$detailed)]) # or select pitch, HNR, ...
feature_mat = t(apply(feature_mat, 1, scale))  # z-transform all variables
feature_mat[is.na(feature_mat)] = 0  # get rid of NAs
colnames(feature_mat) = an$detailed$time  # time stamps in ms
rownames(feature_mat) = 1:nrow(feature_mat)
image(t(feature_mat))  # not a spectrogram, just a feature matrix
getSurprisal(sound, 16000, specFun = feature_mat, logSpec = FALSE)

# analyze all sounds in a folder
surp = getSurprisal('~/Downloads/temp/', savePlots = TRUE)
surp$summary

## End(Not run)

soundgen documentation built on Sept. 20, 2026, 5:07 p.m.