OSCR

Measuring surprisal in sound sequences.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches
  1. [1] § General discussion ↔ analysis/SSM.R, lines 1–106 · score 0.54 · auditory spectrograms, syllable, tune, STFT, segmentation, SSM
  2. [2] § General discussion ↔ analysis/surprisal.R, lines 421–504 · score 0.50 · amplitude envelope, auditory spectrograms, tune, windows, surprisal

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 692 lines · 24 KB · no license · 1 match

  1. #' Self-similarity matrix
  2. #'
  3. #' Calculates the self-similarity matrix and novelty vector of a sound. The
  4. #' self-similarity matrix is produced by cross-correlating different segments of
  5. #' the input sound. Novelty is calculated by convolving the self-similarity
  6. #' matrix with a tapered checkerboard kernel. The positive lobes of the kernel
  7. #' represent coherence (self-similarity within the regions on either side of the
  8. #' center point) and the negative lobes anti-coherence (cross-similarity
  9. #' between these two regions). Since novelty is the dot product of the
  10. #' checkerboard kernel with the SSM, it is high when the two regions are
  11. #' self-similar (internally consistent) but different from each other.
  12. #'
  13. #' @seealso \code{\link{spectrogram}} \code{\link{modulationSpectrum}}
  14. #' \code{\link{segment}}
  15. #'
  16. #' @references \itemize{
  17. #' \item Foote, J. (1999, October). Visualizing music and
  18. #' audio using self-similarity. In Proceedings of the seventh ACM
  19. #' international conference on Multimedia (Part 1) (pp. 77-80). ACM.
  20. #' \item
  21. #' Foote, J. (2000). Automatic audio segmentation using a measure of audio
  22. #' novelty. In Multimedia and Expo, 2000. ICME 2000. 2000 IEEE International
  23. #' Conference on (Vol. 1, pp. 452-455). IEEE.
  24. #' }
  25. #' @inheritParams spectrogram
  26. #' @inheritParams analyze
  27. #' @param sparse if TRUE, the entire SSM is not calculated, but only the central
  28. #' region needed to extract the novelty contour (speeds up the processing)
  29. #' @param input the spectral representation used to calculate the SSM: "audSpec"
  30. #' = auditory spectrogram returned by \code{\link{audSpectrogram}}, "mfcc" =
  31. #' Mel-Frequency Cepstral coefficients, "melspec" = Mel-transformed STFT
  32. #' spectrogram, "spec" = STFT power spectrogram (all three returned by
  33. #' \code{\link[tuneR]{melfcc}}). Any custom spectrogram-like matrix of
  34. #' features (time in columns labeled in s, features in rows) is also accepted
  35. #' (see examples)
  36. #' @param takeLog if TRUE, the input is log-transformed prior to calculating
  37. #' self-similarity
  38. #' @param MFCC which mel-frequency cepstral coefficients to use; defaults to
  39. #' \code{2:13}
  40. #' @param melfcc_pars a list of parameters passed to \code{\link[tuneR]{melfcc}}
  41. #' @param audSpec_pars a list of parameters passed to
  42. #' \code{\link{audSpectrogram}} (if input = 'audSpec')
  43. #' @param norm if TRUE, the spectrum of each STFT frame is normalized
  44. #' @param simil method for comparing frames: "cosine" = cosine similarity, "cor"
  45. #' = Pearson's correlation
  46. #' @param kernelLen length of checkerboard kernel for calculating novelty, ms
  47. #' (larger values favor global, slow vs. local, fast novelty)
  48. #' @param kernelSD SD of checkerboard kernel for calculating novelty
  49. #' @param padWith how to treat edges when calculating novelty: NA = treat sound
  50. #' before and after the recording as unknown, 0 = treat it as silence
  51. #' @param ssmWin window for averaging SSM, frames (has a smoothing effect and
  52. #' speeds up the processing)
  53. #' @param output what to return (drop "ssm" to save memory when analyzing a lot
  54. #' of files)
  55. #' @param plot if TRUE, plots the SSM
  56. #' @param heights relative sizes of the SSM and spectrogram/novelty plot
  57. #' @param main plot title
  58. #' @param specPars graphical parameters passed to \code{filled.contour.mod} and
  59. #' affecting the \code{\link{spectrogram}}
  60. #' @param ssmPars graphical parameters passed to \code{filled.contour.mod} and
  61. #' affecting the plot of SSM
  62. #' @param noveltyPars graphical parameters passed to
  63. #' \code{\link[graphics]{lines}} and affecting the novelty contour
  64. #' @return Returns a list of two components: $ssm contains the self-similarity
  65. #' matrix, and $novelty contains the novelty vector.
  66. #' @export
  67. #' @examples
  68. #' sound = c(soundgen(),
  69. #' soundgen(nSyl = 4, sylLen = 50, pauseLen = 70,
  70. #' formants = NA, pitch = c(500, 330)))
  71. #' # playme(sound)
  72. #' # detailed, local features (captures each syllable)
  73. #' s1 = ssm(sound, samplingRate = 16000, kernelLen = 100,
  74. #' sparse = TRUE) # much faster with 'sparse'
  75. #' # more global features (captures the transition b/w the two sounds)
  76. #' s2 = ssm(sound, samplingRate = 16000, kernelLen = 400, sparse = TRUE)
  77. #'
  78. #' s2$summary
  79. #' s2$novelty # novelty contour
  80. #' \dontrun{
  81. #' ssm(sound, samplingRate = 16000,
  82. #' input = 'mfcc', simil = 'cor', norm = TRUE,
  83. #' ssmWin = 10, # speed up the processing
  84. #' kernelLen = 300, # global features
  85. #' specPars = list(colorTheme = 'seewave'),
  86. #' ssmPars = list(col = rainbow(100)),
  87. #' noveltyPars = list(type = 'l', lty = 3, lwd = 2))
  88. #'
  89. #' # Custom input: produce a nice spectrogram first, then feed it into ssm()
  90. #' sp = spectrogram(sound, 16000, windowLength = c(5, 40), contrast = .3,
  91. #' output = 'processed') # return the modified spectrogram
  92. #' colnames(sp) = as.numeric(colnames(sp)) / 1000 # convert ms to s
  93. #' ssm(sound, 16000, kernelLen = 400, input = sp)
  94. #'
  95. #' # Custom input: use acoustic features returned by analyze()
  96. #' an = analyze(sound, 16000, windowLength = 20, novelty = NULL)
  97. #' input_an = t(an$detailed[, 4:ncol(an$detailed)]) # or select pitch, HNR, ...
  98. #' input_an = t(apply(input_an, 1, scale)) # z-transform all variables
  99. #' input_an[is.na(input_an)] = 0 # get rid of NAs
  100. #' colnames(input_an) = an$detailed$time / 1000 # time stamps in s
  101. #' rownames(input_an) = 1:nrow(input_an)
  102. #' image(t(input_an)) # not a spectrogram, just a feature matrix
  103. #' ssm(sound, 16000, kernelLen = 500, input = input_an, takeLog = FALSE,
  104. #' specPars = list(ylab = 'Feature'))
  105. #' }
  106. ssm = function(
  107. x,
  108. samplingRate = NULL,
  109. from = NULL,
  110. to = NULL,
  111. sparse = FALSE,
  112. input = c('melspec', 'mfcc', 'spec', 'audSpec')[1],
  113. melfcc_pars = list(windowLength = 125, step = 25, nbands = 50),
  114. MFCC = 2:13,
  115. audSpec_pars = list(nFilters = 16, step = 10),
  116. takeLog = FALSE,
  117. norm = FALSE,
  118. simil = c('cosine', 'cor')[1],
  119. kernelLen = 1000,
  120. kernelSD = .5,
  121. padWith = 0,
  122. ssmWin = NULL,
  123. summaryFun = c('mean', 'sd'),
  124. output = c('ssm', 'novelty', 'summary'),
  125. reportEvery = NULL,
  126. cores = 1,
  127. plot = TRUE,
  128. savePlots = NULL,
  129. main = NULL,
  130. heights = c(2, 1),
  131. width = 900,
  132. height = 500,
  133. units = 'px',
  134. res = NA,
  135. specPars = list(
  136. colorTheme = c('bw', 'seewave', 'heat.colors', '...')[2],
  137. xlab = 'Time, s'
  138. ),
  139. ssmPars = list(
  140. colorTheme = c('bw', 'seewave', 'heat.colors', '...')[2],
  141. xlab = 'Time, s',
  142. ylab = 'Time, s'
  143. ),
  144. noveltyPars = list(
  145. type = 'b',
  146. pch = 16,
  147. col = 'black',
  148. lwd = 3
  149. )) {
  150. ## Prepare a list of arguments to pass to .ssm()
  151. myPars = as.list(environment())
  152. # exclude unnecessary args
  153. myPars = myPars[!names(myPars) %in% c(
  154. 'x', 'samplingRate', 'from', 'to', 'savePlots', 'reportEvery', 'cores',
  155. 'summaryFun', 'specPars', 'ssmPars', 'noveltyPars', 'ssmWin',
  156. 'melfcc_pars', 'audSpec_pars')]
  157. myPars$specPars = specPars
  158. myPars$ssmPars = ssmPars
  159. myPars$noveltyPars = noveltyPars
  160. myPars$melfcc_pars = melfcc_pars
  161. myPars$audSpec_pars = audSpec_pars
  162. myPars$win = ssmWin
  164. # analyze
  165. pa = processAudio(
  166. x,
  167. samplingRate = samplingRate,
  168. from = from,
  169. to = to,
  170. funToCall = '.ssm',
  171. myPars = myPars,
  172. reportEvery = reportEvery,
  173. cores = cores,
  174. savePlots = savePlots
  175. )
  177. # htmlPlots
  178. if (!is.null(pa$input$savePlots) && pa$input$n > 1) {
  179. try(htmlPlots(pa$input, savePlots = savePlots, changesAudio = FALSE,
  180. suffix = "ssm", width = paste0(width, units)))
  181. }
  183. # prepare output
  184. if (!is.null(summaryFun) && any(!is.na(summaryFun))) {
  185. temp = vector('list', pa$input$n)
  186. for (i in 1:pa$input$n) {
  187. if (!pa$input$failed[i]) {
  188. temp[[i]] = summarizeAnalyze(
  189. data.frame(novelty = pa$result[[i]]$novelty),
  190. summaryFun = summaryFun,
  191. var_noSummary = NULL)
  192. }
  193. }
  194. idx_failed = which(pa$input$failed)
  195. if (length(idx_failed) > 0) {
  196. idx_ok = which(!pa$input$failed)
  197. if (length(idx_ok) > 0) {
  198. filler = temp[[idx_ok[1]]] [1, ]
  199. filler[1, ] = NA
  200. } else {
  201. stop('Failed to analyze any input')
  202. }
  203. for (i in idx_failed) temp[[i]] = filler
  204. }
  205. mysum_all = cbind(data.frame(file = pa$input$filenames_base),
  206. do.call('rbind', temp))
  207. } else {
  208. mysum_all = NULL
  209. }
  211. if (pa$input$n == 1) {
  212. # unlist
  213. ssm = pa$result[[1]]$ssm
  214. novelty = pa$result[[1]]$novelty
  215. } else {
  216. ssm = lapply(pa$result, function(x) x[['ssm']])
  217. novelty = lapply(pa$result, function(x) x[['novelty']])
  218. }
  219. out = list(ssm = ssm, novelty = novelty, summary = mysum_all)
  220. invisible(out[which(names(out) %in% output)])
  221. }
  224. #' SSM per sound
  225. #'
  226. #' Internal soundgen function.
  227. #' @inheritParams ssm
  228. #' @param audio a list returned by \code{readAudio}
  229. #' @keywords internal
  230. .ssm = function(
  231. audio,
  232. sparse = FALSE,
  233. input = c('melspec', 'mfcc', 'spec', 'audSpec')[1],
  234. melfcc_pars = list(windowLength = 125, step = 25, nbands = 50),
  235. MFCC = 2:13,
  236. audSpec_pars = list(nFilters = 16, step = 10),
  237. takeLog = FALSE,
  238. norm = FALSE,
  239. simil = c('cosine', 'cor')[1],
  240. kernelLen = 100,
  241. kernelSD = .5,
  242. padWith = 0,
  243. win = 1,
  244. output = c('ssm', 'novelty', 'summary'),
  245. plot = TRUE,
  246. main = NULL,
  247. heights = c(2, 1),
  248. width = 900,
  249. height = 500,
  250. units = 'px',
  251. res = NA,
  252. specPars = list(
  253. colorTheme = c('bw', 'seewave', 'heat.colors', '...')[2],
  254. xlab = 'Time, s',
  255. ylab = 'kHz'
  256. ),
  257. ssmPars = list(
  258. colorTheme = c('bw', 'seewave', 'heat.colors', '...')[2],
  259. xlab = 'Time, s',
  260. ylab = 'Time, s'
  261. ),
  262. noveltyPars = list(
  263. type = 'b',
  264. pch = 16,
  265. col = 'black',
  266. lwd = 3
  267. )) {
  268. nyquist = audio$samplingRate / 2
  269. if (is.matrix(input)) {
  270. # custom input to SSM - use as is
  271. target_spec = as.matrix(input)
  272. step = diff(as.numeric(colnames(target_spec))[1:2])
  273. frame_points = round(audio$samplingRate * step)
  274. input = 'custom'
  275. } else {
  276. ## compute mel-filtered spectrum and MFCCs
  277. if (input == 'audSpec') {
  278. # call audSpectrogram()
  279. if (is.null(audSpec_pars$nFilters)) audSpec_pars$nFilters = 32
  280. if (is.null(audSpec_pars$step)) audSpec_pars$step = 20
  281. frame_points = round(audio$samplingRate * audSpec_pars$step / 1000)
  283. target_spec = do.call(.audSpectrogram, c(audSpec_pars, list(
  284. audio = audio,
  285. plot = FALSE
  286. )))$audSpec # cols = time, rows = freq
  287. colnames(target_spec) = as.numeric(colnames(target_spec)) / 1000 # ms to s
  288. } else if (input %in% c('mfcc', 'spec', 'melspec')) {
  289. # call tuneR::melfcc()
  290. if (is.null(melfcc_pars$windowLength)) melfcc_pars$windowLength = 25
  291. if (is.null(melfcc_pars$step)) melfcc_pars$step = 5
  292. frame_points = round(audio$samplingRate * melfcc_pars$step / 1000)
  293. if (!is.numeric(melfcc_pars$windowLength) | melfcc_pars$windowLength <= 0 |
  294. melfcc_pars$windowLength > (audio$duration / 2 * 1000)) {
  295. melfcc_pars$windowLength = min(50, round(audio$duration / 2 * 1000))
  296. warning(paste0(
  297. '"windowLength" must be between 0 and half the sound duration (in ms);
  298. resetting to ', melfcc_pars$windowLength, ' ms')
  299. )
  300. }
  301. if (is.null(melfcc_pars$step))
  302. melfcc_pars$step = melfcc_pars$windowLength / 4
  303. if (is.null(melfcc_pars$nbands)) {
  304. melfcc_pars$nbands = round(100 * melfcc_pars$windowLength / 20)
  305. }
  306. windowLength_points = floor(melfcc_pars$windowLength / 1000 *
  307. audio$samplingRate / 2) * 2
  308. if (is.null(melfcc_pars$maxfreq)) {
  309. melfcc_pars$maxfreq = floor(audio$samplingRate / 2) # Nyquist
  310. }
  312. sound = tuneR::Wave(left = audio$sound, samp.rate = audio$samplingRate, bit = 16)
  313. mel = do.call(tuneR::melfcc, c(
  314. melfcc_pars[which(!names(melfcc_pars) %in% c('windowLength', 'step'))],
  315. list(
  316. samples = sound,
  317. wintime = melfcc_pars$windowLength / 1000,
  318. hoptime = melfcc_pars$step / 1000,
  319. spec_out = TRUE,
  320. numcep = max(MFCC)
  321. )))
  322. if (input == 'mfcc') {
  323. # the first cepstrum presumably makes no sense with amplitude normalization
  324. # (?), and it overestimates the similarity of different frames
  325. target_spec = t(mel$cepstra)[MFCC, ]
  326. target_spec[is.na(target_spec)] = 0 # MFCC are NaN for silent frames
  327. colnames(target_spec) = seq(audio$timeShift, audio$duration,
  328. length.out = ncol(target_spec))
  329. rownames(target_spec) = 1:nrow(target_spec)
  330. } else if (input == 'melspec') {
  331. target_spec = t(mel$aspectrum) # cols = time, rows = freq
  332. colnames(target_spec) = seq(audio$timeShift, audio$duration,
  333. length.out = ncol(target_spec))
  334. rownames(target_spec) = otherToHz(
  335. seq(0, HzToOther(nyquist, "mel"),
  336. length.out = nrow(target_spec)), "mel") / 1000
  337. } else if (input == 'spec') {
  338. target_spec = t(mel$pspectrum) # cols = time, rows = freq
  339. colnames(target_spec) = seq(audio$timeShift, audio$duration,
  340. length.out = ncol(target_spec))
  341. rownames(target_spec) = seq(0, nyquist,
  342. length.out = nrow(target_spec)) / 1000
  343. }
  344. }
  345. }
  346. if (takeLog) {
  347. target_spec = target_spec - min(target_spec, na.rm = TRUE)
  348. target_spec = log(target_spec + min(target_spec[target_spec > 0], na.rm = TRUE))
  349. }
  350. # image(t(target_spec))
  352. ## compute self-similarity matrix
  353. # kernel size in frames, guaranteed to be even
  354. kernelSize = max(4, round(kernelLen * audio$samplingRate / 1000 /
  355. frame_points / 2) * 2)
  356. s = selfsim(
  357. m = target_spec,
  358. norm = norm,
  359. simil = simil,
  360. win = win,
  361. sparse = sparse,
  362. kernelSize = kernelSize
  363. )
  364. # s = zeroOne(s^2) # hist(s)
  365. # image(s)
  367. ## compute novelty
  368. novelty = getNovelty(ssm = s, kernelSize = kernelSize,
  369. kernelSD = kernelSD, padWith = padWith)
  371. ## PLOTTING
  372. if (is.character(audio$savePlots)) {
  373. plot = TRUE
  374. png(filename = paste0(audio$savePlots, audio$filename_noExt, "_ssm.png"),
  375. width = width, height = height, units = units, res = res)
  376. }
  377. if (plot) {
  378. if (is.null(main)) {
  379. if (audio$filename_base == 'sound') {
  380. main = ''
  381. } else {
  382. main = audio$filename_base
  383. }
  384. }
  386. op = par(c('mar', 'xaxt', 'yaxt', 'mfrow')) # save user's original pars
  387. layout(matrix(c(2, 1), nrow = 2, byrow = TRUE), heights = heights)
  388. par(mar = c(5.1, 4.1, 0, 2.1),
  389. xaxt = 's',
  390. yaxt = 's')
  392. # spectrogram
  393. if (input == 'audSpec') {
  394. # spec = zeroOne(t(log(target_spec + 1e-4)))
  395. spec = t(target_spec)
  396. specPars1 = list(
  397. colorTheme = 'seewave',
  398. xlab = 'Time, s',
  399. ylab = 'kHz'
  400. )
  401. specPars1[names(specPars)] = specPars
  402. specPars1$color.palette = switchColorTheme(specPars1$colorTheme)
  403. specPars1[['colorTheme']] = NULL
  404. do.call(filled.contour.mod, c(list(
  405. x = as.numeric(rownames(spec)),
  406. y = as.numeric(colnames(spec)),
  407. z = spec,
  408. yScale = if (is.null(audSpec_pars$yScale)) 'bark' else audSpec_pars$yScale
  409. ), specPars1))
  410. specPars1$xlim = range(as.numeric(rownames(spec)))
  411. specPars1$ylim = range(as.numeric(colnames(spec)))
  412. } else {
  413. # # log-transform and normalize spectrogram
  414. # if (input == 'melspec') {
  415. # spec = log(zeroOne(mel$aspectrum) + 1e-4) # dynamic range ~ 80 dB or 1e-4
  416. # } else {
  417. # spec = log(zeroOne(mel$pspectrum) + 1e-4)
  418. # }
  419. # spec = zeroOne(spec)
  420. spec = t(target_spec)
  421. specPars1 = list(
  422. colorTheme = 'seewave',
  423. xlab = 'Time, s',
  424. ylab = if (input %in% c('custom', 'mfcc')) input else 'kHz'
  425. )
  426. specPars1[names(specPars)] = specPars
  427. specPars1$color.palette = switchColorTheme(specPars1$colorTheme)
  428. specPars1[['colorTheme']] = NULL
  429. do.call(filled.contour.mod, c(list(
  430. x = as.numeric(rownames(spec)),
  431. y = as.numeric(colnames(spec)),
  432. z = spec,
  433. yScale = if (input == 'melspec') 'mel' else 'linear'
  434. ), specPars1
  435. ))
  436. }
  437. specPars1$xlim = c(audio$timeShift, audio$duration)
  438. specPars1$ylim = range(as.numeric(rownames(target_spec)))
  439. if (input == 'melspec')
  440. specPars1$ylim = HzToOther(specPars1$ylim * 1000, 'mel')
  442. # novelty
  443. noveltyPars1 = list(
  444. type = 'b',
  445. pch = 16,
  446. col = 'black',
  447. lwd = 3
  448. )
  449. noveltyPars1[names(noveltyPars)] = noveltyPars
  450. do.call(lines, c(list(
  451. x = seq(specPars1$xlim[1], specPars1$xlim[2], length.out = length(novelty)),
  452. y = zeroOne(novelty) * specPars1$ylim[2] * .95
  453. ), noveltyPars1
  454. ))
  455. axis(side = 1, labels = TRUE)
  456. par(mar = c(0, 4.1, 2.1, 2.1),
  457. xaxt = 'n',
  458. yaxt = 's')
  459. xlab = ''
  461. # SSM
  462. ssmPars1 = list(
  463. levels = seq(0, 1, length = 30),
  464. colorTheme = 'seewave',
  465. xlab = 'Time, s',
  466. ylab = 'Time, s',
  467. main = main
  468. )
  469. ssmPars1[names(ssmPars)] = ssmPars
  470. ssmPars1$color.palette = switchColorTheme(ssmPars1$colorTheme)
  471. ssmPars1[['colorTheme']] = NULL
  472. timestamps_ssm = seq(0, audio$duration, length.out = nrow(s))
  473. do.call(filled.contour.mod, c(list(
  474. x = timestamps_ssm,
  475. y = timestamps_ssm,
  476. z = s,
  477. y_Hz = FALSE
  478. ), ssmPars1
  479. ))
  480. # restore original pars
  481. par('mar' = op$mar, 'xaxt' = op$xaxt, 'yaxt' = op$yaxt, 'mfrow' = op$mfrow)
  482. if (is.character(audio$savePlots)) dev.off()
  483. }
  485. out = list(ssm = s, novelty = novelty)
  486. invisible(out[which(names(out) %in% output)])
  487. }
  490. #' Compute self-similarity
  491. #'
  492. #' Internal soundgen function.
  493. #'
  494. #' Called by \code{\link{ssm}}.
  495. #' @param m input matrix such as a spectrogram
  496. #' @inheritParams ssm
  497. #' @param win the length of window for averaging self-similarity, frames
  498. #' @return Returns a square self-similarity matrix.
  499. #' @keywords internal
  500. #' @examples
  501. #' m = matrix(rnorm(40), nrow = 5)
  502. #' soundgen:::selfsim(m, sparse = TRUE, kernelSize = 2)
  503. selfsim = function(m,
  504. norm = FALSE,
  505. simil = c('cosine', 'cor')[1],
  506. win = 1,
  507. sparse = FALSE,
  508. kernelSize = NULL) {
  509. nc = ncol(m)
  510. if (win > floor(nc / 2)) {
  511. win = floor(nc / 2)
  512. warning(paste('"win" must be smaller than half the number of frames',
  513. 'resetting to', floor(nc / 2)))
  514. }
  515. if (win %% 2 == 0) {
  516. win = max(ceiling(win / 2) * 2 - 1, 1)
  517. } # win must be odd
  519. # normalize input by column, if needed
  520. if (norm) {
  521. m = apply(m, 2, zeroOne, na.rm = TRUE)
  522. m[is.na(m)] = 0
  523. }
  525. # calculate windows for averaging self-similarity
  526. winIdx = unique(round(seq(1, nc - win + 1, length.out = ceiling(nc / win))))
  527. numWins = length(winIdx)
  529. # calculate the lower triangle of self-similarity matrix
  530. out = matrix(NA, nrow = numWins, ncol = numWins)
  531. rownames(out) = colnames(out) = winIdx
  532. if (!sparse) j_idx = seq_len(numWins)
  533. for (i in seq_along(winIdx)) {
  534. if (sparse) {
  535. j_idx = max(1, i - kernelSize) : max(1, (i - 1))
  536. } else {
  537. j_idx = 1:max(1, (i - 1))
  538. }
  539. for (j in j_idx) {
  540. mi = as.vector(m[, winIdx[i]:(winIdx[i] + win - 1)])
  541. mj = as.vector(m[, winIdx[j]:(winIdx[j] + win - 1)])
  542. if (any(mi != 0) && any(mj != 0)) {
  543. if (simil == 'cosine') {
  544. # http://stackoverflow.com/questions/6597005/cosine-similarity-between-two-vectors-in-language-r
  545. out[i, j] = crossprod(mi, mj) / sqrt(crossprod(mi) * crossprod(mj))
  546. } else if (simil == 'cor') {
  547. out[i, j] = cor(mi, mj)
  548. }
  549. } else {
  550. # if at least one is a vector of zeros, set result to 0 (otherwise NA)
  551. out[i, j] = 0
  552. }
  553. }
  554. }
  555. # fill up the upper triangle as well
  556. diag(out) = 1
  557. out1 = t(out)
  558. out1[lower.tri(out1)] = out[lower.tri(out)]
  559. # isSymmetric(out1)
  560. # image(out1)
  561. zeroOne(t(out1), na.rm = TRUE)
  562. }
  565. #' Checkerboard kernel
  566. #'
  567. #' Internal soundgen function.
  568. #'
  569. #' Prepares a square matrix \code{size x size} specifying a gaussian kernel for
  570. #' measuring novelty of self-similarity matrices. Called by
  571. #' \code{\link{getNovelty}}
  572. #' @param size kernel size (points), preferably an even number
  573. #' @param kernel_mean,kernelSD mean and SD of the gaussian kernel
  574. #' @param plot if TRUE, shows a perspective plot of the kernel
  575. #' @param checker if TRUE, inverts two quadrants
  576. #' @return Returns a square matrix with \code{size} rows and columns.
  577. #' @keywords internal
  578. #' @examples
  579. #' kernel = soundgen:::getCheckerboardKernel(size = 64, kernelSD = 0.1, plot = TRUE)
  580. #' dim(kernel)
  581. #' kernel = soundgen:::getCheckerboardKernel(size = 19, kernelSD = .5,
  582. #' checker = FALSE, plot = TRUE)
  583. #' kernel = soundgen:::getCheckerboardKernel(size = c(9, 45), kernelSD = .5,
  584. #' checker = FALSE, plot = TRUE)
  585. #' kernel = soundgen:::getCheckerboardKernel(size = c(9, 45), kernelSD = .5,
  586. #' checker = TRUE, plot = TRUE)
  587. getCheckerboardKernel = function(size,
  588. kernel_mean = 0,
  589. kernelSD = 0.5,
  590. plot = FALSE,
  591. checker = TRUE) {
  592. if (length(size) == 1) {
  593. x = y = seq(-1, 1, length.out = size)
  594. size = c(size, size)
  595. } else if (length(size) == 2) {
  596. x = seq(-1, 1, length.out = size[1])
  597. y = seq(-1, 1, length.out = size[2])
  598. } else {
  599. stop('size must be of length 1 or 2')
  600. }
  602. kernelSD = kernelSD # just to get rid of the "unused arg" warning in CMD check :-)
  603. if (max(size) < 50) {
  604. # faster than mvtnorm::dmvnorm for small kernels
  605. kernel = matrix(NA, ncol = size[2], nrow = size[1])
  606. for (i in seq_len(nrow(kernel))) {
  607. for (j in seq_len(ncol(kernel))) {
  608. kernel[i, j] = dnorm(x[i], mean = kernel_mean, sd = kernelSD) *
  609. dnorm(y[j], mean = kernel_mean, sd = kernelSD)
  610. }
  611. }
  612. } else {
  613. # this is faster for large kernels
  614. sigma = diag(2) * kernelSD
  615. kernel_long = expand.grid(x1 = x, x2 = y)
  616. kernel_long$dd = mvtnorm::dmvnorm(x = kernel_long,
  617. mean = c(kernel_mean, kernel_mean),
  618. sigma = sigma)
  619. kernel = matrix(kernel_long$dd, nrow = size[1])
  620. # kernel[1:5, 1:5]
  621. }
  623. if (checker) {
  624. fl_row = floor(size[1] / 2)
  625. fl_col = floor(size[2] / 2)
  626. cl_row = ceiling(size[1] / 2)
  627. cl_col = ceiling(size[2] / 2)
  628. # quadrant 0 to 3 o'clock
  629. kernel[seq_len(fl_row), (cl_col + 1):size[2]] = -kernel[seq_len(fl_row), (cl_col + 1):size[2]]
  630. # quadrant 6 to 9 o'clock
  631. kernel[(cl_row + 1):size[1], seq_len(cl_col)] = -kernel[(cl_row + 1):size[1], seq_len(cl_col)]
  632. }
  633. kernel = kernel / max(kernel)
  635. if (plot) {
  636. persp(
  637. kernel,
  638. theta = -20,
  639. phi = 25,
  640. # zlim = c(-1, 4),
  641. ticktype = 'detailed'
  642. )
  643. }
  644. kernel
  645. }
  648. #' SSM novelty
  649. #'
  650. #' Internal soundgen function.
  651. #'
  652. #' Calculates novelty in a self-similarity matrix. Called by \code{\link{ssm}}.
  653. #' @param ssm self-similarity matrix, as produced by \code{\link{selfsim}}
  654. #' @param kernelSize the size of gausisan kernel (points)
  655. #' @param kernelSD the SD of gaussian kernel
  656. #' @param normalize if TRUE, normalizes so that max = 1
  657. #' @return Returns a numeric vector of length \code{nrow(ssm)}
  658. #' @keywords internal
  659. getNovelty = function(ssm,
  660. kernelSize,
  661. kernelSD,
  662. padWith = 0,
  663. normalize = TRUE) {
  664. kernel = getCheckerboardKernel(size = kernelSize, kernelSD = kernelSD)
  665. ## pad matrix with size / 2 zeros, so that we can correlate it with the
  666. # kernel starting from the very edge
  667. ssm_padded = matrix(padWith,
  668. nrow = nrow(ssm) + kernelSize,
  669. ncol = nrow(ssm) + kernelSize)
  670. halfK = kernelSize / 2
  671. # indices in the padded matrix where we'll paste the original ssm
  672. idx = c(halfK + 1, nrow(ssm_padded) - halfK)
  673. # paste original. Now we have a padded ssm
  674. ssm_padded[idx[1]:idx[2], idx[1]:idx[2]] = ssm
  676. ## get novelty
  677. novelty = rep(0, nrow(ssm))
  678. # for each point on the main diagonal, novelty = correlation between the
  679. # checkerboard kernel and the ssm. See Badawy, "Audio novelty-based
  680. # segmentation of music concerts"
  681. for (i in idx[1]:idx[2]) {
  682. n = (i - halfK):(i + halfK - 1)
  683. # suppress warnings, b/c otherwise cor complains of sd = 0 for silent segments
  684. mat_i = ssm_padded[n, n]
  685. diag(mat_i) = NA
  686. novelty[i - halfK] = suppressWarnings(
  687. cor(as.vector(mat_i), as.vector(kernel), use = 'pairwise.complete.obs')
  688. )
  689. }
  690. novelty[is.na(novelty)] = 0
  691. novelty
  692. }

SSM.R, no license · at the source

Overview

Authors: Andrey Anikin1
ORCID iDs: Andrey Anikin
  1. Division of Cognitive Science, Department of Philosophy, Lund University, Box 192, SE-221 00 Lund, Sweden
Institutions: Lund University (Sweden)
Journal: Behavior research methods, volume 58, issue 10, article 278
Dates: received 13 February 2026; accepted 21 July 2026; published online 24 August 2026; in print 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.3758/s13428-026-03153-3 · PMID 42637973 · PMCID PMC13503466 · OpenAlex W7128612809
Open access: hybrid, a free copy (OpenAlex)
Preprint: osf.io/bgzvc
Status: code verified
Categories: human (organism)
Keywords: Auditory attention, Salience, Shannon surprisal, Bayesian surprise, Self-similarity
MeSH: Acoustic Stimulation*, Auditory Perception*, Sound*, Algorithms, Animals, Attention, Bayes Theorem, Humans, Neural Networks, Computer (* major topic)
Topic: Hearing Loss and Rehabilitation (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Vetenskapsrådet (2023-00850)
Citations: not cited yet (Europe PMC); 68 references in the paper

Abstract

Sensory input that violates prior expectations attracts attention, making unpredictability an important perceptual property to measure. In the auditory modality, knowing what sounds will be perceived as surprising, and therefore salient, is relevant both for studying vocal communication and for applied purposes such as managing noise pollution. Focusing on sequences of animal vocalizations and environmental sounds as ecologically important acoustic stimuli, I describe and benchmark several algorithms for measuring their perceived unpredictability. Information-theoretical approaches include Shannon surprisal and Bayesian surprise, both implemented here to detect deviant stimuli based on distributional acoustic properties. The second group of algorithms is based on detecting spectro-temporal recurrence assessed with autocorrelation functions (ACF surprisal) and self-similarity matrices (SSM novelty). The third approach uses neural networks. Based on the ratings of the predictability of 300 synthetic acoustic sequences by 195 human listeners, Shannon surprisal and SSM novelty capture the perceived unpredictability that is due to spectral variability, whereas ACF surprisal taps into the perceptual impact of irregular rhythm. Most algorithms converge on the time scale of about 1 s as the most perceptually relevant for spectral variability, which is consistent with the hypothesis that the perception of unpredictability stems from a relatively limited amount of auditory input held in short-term memory. Together, the presented open-source algorithms offer powerful and flexible tools for measuring acoustic surprisal and studying auditory attention, while the corpus of predictability ratings offers a resource for future benchmarking. All code and data are freely available from the R package soundgen and supplementary materials at https://osf.io/bgzvc.

Supplementary Information: The online version contains supplementary material available at 10.3758/s13428-026-03153-3.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

OSF bgzvc

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Languages: R (16)
Size: 54 files, 16 scripts
Software Heritage: not checked
Found in: the text, “Data analysis”
Holds: 6 notebooks
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: ggplot2 (7 files), brms (5 files), patchwork (1 file), randomForest (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
16 files
At the source: osf.io/bgzvc

OSF kp9mg

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Languages: R (16)
Size: 42 files, 16 scripts
Software Heritage: not checked
Found in: the references
Holds: 5 notebooks
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: brms (3 files), ggplot2 (3 files), tidyverse (3 files), reshape2 (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
16 files
At the source: osf.io/kp9mg/

Code availability

All data and materials are available at https://osf.io/bgzvc. This includes R code for running all reported analyses. The functions for measuring surprisal are included in the open-source R package soundgen, freely available from https://cran.r-project.org/package=soundgen.

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 32 scripts, each with its path and the digest of its content;
  • 2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability Statement

All data and materials are available at https://osf.io/bgzvc.

All data and materials are available at https://osf.io/bgzvc. This includes R code for running all reported analyses. The functions for measuring surprisal are included in the open-source R package soundgen, freely available from https://cran.r-project.org/package=soundgen.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

  • Publisher: n/a → Springer Science+Business Media

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 1 author, 5 keywords, 9 MeSH terms, 1 funder, 44 references.

Cite

This paper

Anikin, A. (2026). Measuring surprisal in sound sequences. Behavior research methods, 58(10), 278. https://doi.org/10.3758/s13428-026-03153-3

BibTeX

@article{anikin2026measuring,
author = {Anikin, Andrey},
title = {{Measuring surprisal in sound sequences}},
journal = {Behavior research methods},
year = {2026},
month = aug,
volume = {58},
number = {10},
pages = {278},
publisher = {Springer Science+Business Media},
issn = {1554-351X},
doi = {10.3758/s13428-026-03153-3},
url = {https://doi.org/10.3758/s13428-026-03153-3},
pmid = {42637973},
pmcid = {PMC13503466}
}

RIS

TY - JOUR
AU - Anikin, Andrey
TI - Measuring surprisal in sound sequences
T2 - Behavior research methods
J2 - Behav Res Methods
PY - 2026
DA - 2026/08/24
VL - 58
IS - 10
SP - 278
SN - 1554-351X
PB - Springer Science+Business Media
DO - 10.3758/s13428-026-03153-3
UR - https://doi.org/10.3758/s13428-026-03153-3
LA - en
ER -

CSL-JSON

{
"id": "10.3758/s13428-026-03153-3",
"type": "article-journal",
"title": "Measuring surprisal in sound sequences",
"container-title": "Behavior research methods",
"author": [
{
"family": "Anikin",
"given": "Andrey"
}
],
"container-title-short": "Behav Res Methods",
"volume": "58",
"issue": "10",
"page": "278",
"DOI": "10.3758/s13428-026-03153-3",
"PMID": "42637973",
"PMCID": "PMC13503466",
"ISSN": "1554-351X",
"publisher": "Springer Science+Business Media",
"URL": "https://doi.org/10.3758/s13428-026-03153-3",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
24
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: randomForest, brms, reshape2, 3 other tools
[2] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: randomForest, reshape2, patchwork, 2 other tools
[3] doi:10.3390/ijms27156925 [code]
XGBoost-SHAP Interpretable Modeling Identifies and Validates an Eight-Gene Biomarker for Hepatic Encephalopathy Risk Prediction in Cirrhosis.
Journal: International journal of molecular sciences
In common: randomForest, reshape2, patchwork, 2 other tools
[4] doi:10.1162/imag.a.1347 [code]
Neural and behavioural correlates of theory of mind reasoning in five-year-old children born preterm.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: brms, reshape2, patchwork, 2 other tools
[5] doi:10.1162/opmi.a.372 [code]
Broadening the Agent Preference Hypothesis Through Experiencers: Eye-Tracking and EEG Evidence of Proto-Agents and Proto-Patients.
Journal: Open mind : discoveries in cognitive science
In common: brms, reshape2, patchwork, 2 other tools
[6] doi:10.1093/braincomms/fcag217 [code]
Investigation of stress hormones across multiday seizure cycles.
Journal: Brain communications
In common: brms, reshape2, patchwork, 2 other tools
[7] doi:10.1038/s41593-026-02387-w [code]
Developing mouse inhibitory neuron single-cell transcriptomes reveal distinct modes of cell-type diversification.
Journal: Nature neuroscience
In common: randomForest, reshape2, patchwork, 2 other tools
[8] doi:10.1016/j.celrep.2026.117505 [code]
Impaired spatial coding and neuronal hyperactivity in the medial entorhinal cortex of aged APP knock-in mice.
Journal: Cell reports
In common: brms, reshape2, patchwork, 2 other tools
[9] doi:10.1162/imag.a.1258 [code]
Non-specific increase in alpha power during a neurofeedback session targeting its downregulation.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: brms, reshape2, patchwork, 2 other tools
[10] doi:10.1038/s44271-026-00431-w [code]
Alpha power increases spontaneously during a neurofeedback session.
Journal: Communications psychology
In common: brms, reshape2, patchwork, 2 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.