Package {soundgen}

The soundgen package

Description

Parametric source-filter synthesis of harmonic-noise signals, such as animal vocalizations and human voice, with control over pitch, formants, noise, amplitude modulation, nonlinear phenomena, and morphing. General signal processing tools for audio analysis and manipulation: pitch tracking, formant and vocal tract length estimation, reassigned and auditory spectrograms, modulation spectra and psychoacoustic roughness, self-similarity and surprisal, audio segmentation, pitch and formant shifting, etc. Includes four interactive web apps for audio synthesis, annotation, formant analysis, and manually correcting pitch contours.

Details

Visualization

• spectrogram Ordinary (STFT) spectrograms

• audSpectrogram Auditory spectrograms

• modulationSpectrum Modulation spectra and roughness (see Anikin 2025)

• ssm Self-similarity matrices, acoustic novelty (see Anikin 2026)

• phasegram Nonlinear analysis, phasegrams (see Anikin & Herbst 2025)

Acoustic analysis

• analyze General-purpose acoustic analysis (extracts nearly everything: spectral descriptives, pitch, formants, loudness, roughness, ...)

• segment Audio segmentation, analysis of rhythm and temporal structure of sequences

• getSurprisal Shannon and autocorrelation-based surprisal, Bayesian surprise (see Anikin 2026)

• estimateVTL, schwa Formant analysis (see Anikin et al. 2024)

Interactive acoustic analysis

• pitch_app Manually corrected formant measurement

• formant_app Manually corrected pitch tracking

• annotation_app Audio annotation

Sound synthesis

• soundgen The main function for sound and voice synthesis (see Anikin 2019)

• soundgen_app Interactive version of the soundgen() function

• morph Morphing of two sounds based on their soundgen code

Sound modification

• normalizeFolder Peak / RMS etc normalization of recordings

• timeStretch, shiftPitch, prosody, shiftFormants Selective modification of fundamental frequency, duration, and/or resonance frequency

• noiseRemoval, compressor, addAM, resample, etc: various audio effects

Vignettes

• https://cogsci.se/soundgen/spectrograms.html

• https://cogsci.se/soundgen/sound_generation.html

• https://cogsci.se/soundgen/acoustic_analysis.html

• https://cogsci.se/soundgen/matching/matching.html

Author(s)

Maintainer: Andrey Anikin andrey.anikin@cogsci.se

Authors:

• Andrey Anikin andrey.anikin@cogsci.se

References

• Anikin, A. (2019). Soundgen: an open-source tool for synthesizing nonverbal vocalizations. Behavior Research Methods, 51(2), 778-792.

• Anikin, A., Barreda, S. & Reby, D. (2024) A practical guide to calculating vocal tract length and scale-invariant formant patterns. Behavior Research Methods 56, 5588–5604.

• Anikin, A. (2025) Acoustic estimation of voice roughness. Attention, Perception, & Psychophysics 87: 1771–1787.

• Anikin, A. & Herbst, C.T. (2025) How to analyze and manipulate nonlinear phenomena in voice recordings. Philosophical Transactions B 380: 20240003.

• Anikin. A. (2026) Measuring surprisal in sound sequences. Behavior Research Methods.

See Also

Useful links:

• http://cogsci.se/soundgen.html

Canonical roxygen docs (inherited, not user-facing)

Description

Canonical roxygen docs (inherited, not user-facing)

Usage

.roxygen_defaults(
  x,
  samplingRate,
  scale,
  from,
  to,
  windowLength,
  step,
  overlap,
  wn,
  zp,
  dynamicRange,
  play,
  saveAudio,
  plot,
  savePlots,
  embed,
  reportEvery,
  cores,
  width,
  height,
  units,
  res,
  normalize
)

Arguments

Convert between Hz and other frequency scales

Description

Converts between Hz and ERB, bark, mel, log, semitones relative to a reference frequency, or musical notes. Accepts vectors and missing values.

Usage

HzToOther(
  x,
  scale = c("ERB", "bark", "mel", "log", "semitones", "notes", "linear", "orig"),
  ...
)

otherToHz(
  x,
  scale = c("ERB", "bark", "mel", "log", "semitones", "notes", "linear", "orig"),
  ...
)

Arguments

Value

Numeric vector of converted frequencies.

Scale‑specific arguments

The following arguments can be passed via ... . They are only used when the corresponding scale is selected; otherwise, they are silently ignored.

method

Character, used when scale = "ERB". Choose between "linear" (Glasberg & Moore, 1990) and "quadratic" (Moore & Glasberg, 1983). Default is "linear".

ref

Numeric, used when scale = "semitones". Reference frequency (Hz) for the semitone scale. Default is C-5 (0.5109875 Hz).

A4

Numeric, used when scale = "notes". Frequency of the note A4. Default is 440 Hz (modern standard ISO 16, concert pitch).

showCents

Logical, used only with HzToOther() when scale = "notes". If TRUE, appends the cent deviation (rounded to the nearest integer) to the note name. Default is FALSE. Not applicable to otherToHz().

No additional arguments are required or usedor scale = "bark", "mel", "log", or "linear".

References

• Moore, B. C., & Glasberg, B. R. (1983). Suggested formulae for calculating auditory-filter bandwidths and excitation patterns. The journal of the acoustical society of America, 74(3), 750-753.

• Glasberg, B. R., & Moore, B. C. (1990). Derivation of auditory filter shapes from notched-noise data. Hearing research, 47(1-2), 103-138.

Examples

x = c(-20, 20, 100, 440, 1000, NA)
HzToOther(x, 'ERB')
HzToOther(x, 'ERB', 'quadratic')
HzToOther(x, 'bark')
HzToOther(x, 'mel')
HzToOther(x, 'log')
HzToOther(x, 'semitones', ref = 16)
HzToOther(x, 'notes', showCents = TRUE)

# ...and back to Hz
x = c(0:10, NA)

otherToHz(x, 'ERB')
otherToHz(x, 'ERB', method = 'quadratic')
otherToHz(HzToOther(c(100, 440, 2000), 'ERB'), 'ERB')

otherToHz(x, 'bark')
otherToHz(HzToOther(c(100, 440, 2000), 'bark'), 'bark')

otherToHz(x, 'mel')
otherToHz(HzToOther(c(100, 440, 2000), 'mel'), 'mel')

otherToHz(x, 'log')

otherToHz(x, 'semitones')
HzToOther(c(440, 210, 880), 'semitones', ref = 440)
otherToHz(HzToOther(c(440, 210, 880), 'semitones'), 'semitones')

otherToHz(c('A4', 'C#6', 'blabla', 'C0', 'C-2', NA), 'notes')
HzToOther(c(440, 293, 115, 16.35, 4), 'notes')

HzToOther(c(440, 415, 80, 81), 'notes', showCents = TRUE)
# 80 Hz is almost exactly midway (+49 cents) between D#2 and E2

# Baroque tuning A415, half a semitone flat relative to concert pitch A440
HzToOther(c(440, 415, 16.35), 'notes', A4 = 415)
otherToHz(c("A4", "D4", "A#2", "C0", "C-2"), 'notes', A4 = 415)

Impulse response

Description

A finger snap in a room as an approximation to an impulse response

Usage

IR

Format

An object of class Wave of length 3252.

Autocorrelation with FFT

Description

Analogous to acf, but based on the Fast Fourier Transform (Wiener-Khinchin theorem) and 5-10 times faster, especially for long input vectors.

Usage

acf_fft(x, center = TRUE, padToMult = 2)

Arguments

Value

Numeric vector that is usually longer than input because it is padded with zeros to the next power of two. Constant and zero inputs return NA.

Examples

len = 200
x = sin(2 * pi * 100 * (1:len) / 1000) + rnorm(len, 0, .5)
plot(x, type = 'l')
aut = acf(x, lag.max = len/2)
aut2 = acf_fft(x)
points(0:100, aut2[1:(len/2+1)], type = 'l', col = 'blue')
aut$acf[1:10]
aut2[1:10]

# compare execution time
system.time(for (i in 1:100) acf(x, lag.max = len/2, plot = FALSE))
system.time(for (i in 1:100) acf_fft(x))

Add amplitude modulation

Description

Adds sinusoidal or logistic amplitude modulation to a sound. Sinusoidal AM creates a single pair of sidebands at ±amFreq around each original harmonic, whereas non-sinusoidal AM creates broader sidebands with more extra harmonics (see examples).

Usage

addAM(
  x,
  samplingRate = NULL,
  amDep = 25,
  amFreq = 30,
  amType = c("logistic", "sine"),
  amShape = 0,
  invalidArgAction = c("adjust", "abort", "ignore"),
  play = FALSE,
  saveAudio = FALSE,
  plot = FALSE,
  savePlots = FALSE,
  embed = FALSE,
  reportEvery = NULL,
  cores = 1,
  width = 900,
  height = 500,
  units = "px",
  res = NA
)

Arguments

Value

Returns the modified audio as a numeric vector with the original sampling rate.

Examples

sound1 = soundgen(pitch = c(200, 300), addSilence = 0)
s1 = addAM(sound1, 16000, amDep = c(0, 50, 0), amFreq = 75, plot = TRUE)
# playme(s1)
## Not run: 
# Parameters can be specified as in the soundgen() function, eg:
s2 = addAM(sound1, 16000,
         amDep = list(time = c(0, 50, 52, 200, 201, 300),
                      value = c(0, 0, 35, 25, 0, 0)),
         plot = TRUE, play = TRUE)

# Sinusoidal AM produces exactly 2 extra harmonics at ±amFreq
# around each f0 harmonic (amFreq, and thus the width of sidebands,
# may vary over time):
s3 = addAM(sound1, 16000, amDep = 30, amFreq = c(50, 80),
           amType = 'sine', plot = TRUE, play = TRUE)
spectrogram(s3, 16000, windowLength = 150, ylim = c(0, 2))

# Non-sinusoidal AM produces multiple new harmonics,
# which can resemble subharmonics...
s4 = addAM(sound1, 16000, amDep = 70, amFreq = 50, amShape = -1,
           plot = TRUE, play = TRUE)
spectrogram(s4, 16000, windowLength = 150, ylim = c(0, 2))

# ...but more often look like sidebands
sound3 = soundgen(sylLen = 600, pitch = c(800, 1300, 1100), addSilence = 0)
s5 = addAM(sound3, 16000, amDep = c(0, 30, 100, 40, 0),
           amFreq = 105, amShape = -.3,
           plot = TRUE, play = TRUE)
spectrogram(s5, 16000, ylim = c(0, 5))

# Feel free to add AM stochastically:
s6 = addAM(sound1, 16000,
           amDep = rnorm(10, 40, 20), amFreq = rnorm(20, 70, 20),
           plot = TRUE, play = TRUE)
spectrogram(s6, 16000, windowLength = 150, ylim = c(0, 2))

# If amFreq is locked to an integer ratio of f0, we can get subharmonics
# For ex., here is with pitch 400-600-400 Hz (soundgen interpolates pitch
# on a log scale and amFreq on a linear scale, so we align them by extracting
# a long contour on a log scale for both)
con = soundgen:::getSmoothContour(anchors = c(400, 600, 400),
                       len = 20, thisIsPitch = TRUE)
s = soundgen(sylLen = 1500, pitch = con, amFreq = con/3, amDep = 30,
             plot = TRUE, play = TRUE, ylim = c(0, 3))

# Process all files in a folder and save the modified audio
addAM('~/Downloads/temp', saveAudio = TRUE, amFreq = 70, amDep = c(0, 50))

## End(Not run)

Add formants

Description

A spectral filter that either adds or removes formants from a sound - that is, amplifies or dampens certain frequency bands, as in human vowels. See soundgen and getFormantFilter for more information. With action = 'remove' this function can perform inverse filtering to remove formants and obtain raw glottal output, provided that you can specify the correct formant structure. Instead of formants, any arbitrary spectral filtering function can be applied using the formantFilter argument (e.g., for a low/high/bandpass filter).

Usage

addFormants(
  x,
  samplingRate = NULL,
  formants = NULL,
  formantFilter = NULL,
  action = c("add", "remove"),
  dB = NULL,
  specificity = 1,
  zFun = NULL,
  vocalTract = NA,
  formantDep = 1,
  formantDepStoch = 1,
  formantWidth = 1,
  formantCeiling = NULL,
  lipRad = 6,
  noseRad = 4,
  mouthOpenThres = 0,
  mouth = NA,
  temperature = 0.025,
  formDrift = 0.3,
  formDisp = 0.2,
  smoothing = list(interpol = "splineFC"),
  windowLength = 50,
  step = NULL,
  overlap = 75,
  wn = "gaussian",
  normalize = c("orig", "max", "none"),
  play = FALSE,
  saveAudio = FALSE,
  reportEvery = NULL,
  cores = 1,
  ...
)

Arguments

Details

Algorithm: converts input from a time series (time domain) to a spectrogram (frequency domain) through short-time Fourier transform (STFT), multiplies by the spectral filter containing the specified formants, and transforms back to a time series via inverse STFT. This is a subroutine for voice synthesis in soundgen, but it can also be applied to a recording.

Value

The filtered waveform as a numeric vector of the original length with the original sampling rate, or a list if there are multiple inputs.

See Also

getFormantFilter transplantFormants soundgen

Examples

sound = c(rep(0, 1000), rnorm(8000) * 2 - 1, rep(0, 1000))  # white noise
# NB: pad with silence to avoid artifacts if removing formants
# playme(sound)
# spectrogram(sound, samplingRate = 16000)

# add F1 = 900, F2 = 1300 Hz
sound_filtered = addFormants(sound, samplingRate = 16000,
                             formants = c(900, 1300))
# playme(sound_filtered)
# spectrogram(sound_filtered, samplingRate = 16000)

# ...and remove them again (assuming we know what the formants are)
sound_inverse_filt = addFormants(sound_filtered,
                                 samplingRate = 16000,
                                 formants = c(900, 1300),
                                 action = 'remove')
# playme(sound_inverse_filt)
# spectrogram(sound_inverse_filt, samplingRate = 16000)

## Not run: 
## Perform some user-defined manipulation of the spectrogram with zFun
# Ex.: noise removal - silence all bins 50 dB below the max value
s_noisy = soundgen(sylLen = 200, addSilence = 0,
                   noise = list(time = c(-100, 300), value = -20))
spectrogram(s_noisy, 16000)
# playme(s_noisy)
zFun = function(z, cutoff = -50) {
  az = abs(z)
  thres = max(az) * 10 ^ (cutoff / 20)
  z[which(az < thres)] = 0
  return(z)
}
s_denoised = addFormants(s_noisy, samplingRate = 16000,
                         formants = NA, zFun = zFun, cutoff = -40)
spectrogram(s_denoised, 16000)
# playme(s_denoised)

# If neither formants nor formantFilter are defined, only lipRad has an effect
# For ex., we can boost low frequencies by 6 dB/oct
noise = rnorm(8000)
noise1 = addFormants(noise, 16000, lipRad = -6)
meanSpectrum(noise1, 16000, yScale = 'max0')

# Arbitrary spectra can be defined with formantFilter. For ex., we can
# have a flat spectrum up to 2 kHz (Nyquist / 4) and -3 dB/kHz above:
freqs = seq(0, 16000 / 2, length.out = 100)
n = length(freqs)
idx = (n / 4):n
sp_dB = c(rep(0, n / 4 - 1), (freqs[idx] - freqs[idx[1]]) / 1000 * (-3))
plot(freqs, sp_dB, type = 'b')
noise2 = addFormants(noise, 16000, lipRad = 0, formantFilter = 10 ^ (sp_dB / 20))
meanSpectrum(noise2, 16000, yScale = 'max0')

## Use the spectral envelope of another recording
# (NB: this can also be achieved with a single call to transplantFormants)
sound_orig = soundgen(sylLen = 300, formants = 'a', addSilence = 5)
samplingRate = 16000
# playme(sound_orig, samplingRate)

# get a few pitch anchors to reproduce the original intonation
pitch = analyze(sound_orig, samplingRate = samplingRate,
  pitchMethod = c('autocor', 'dom'))$detailed$pitch
pitch = pitch[!is.na(pitch)]

# extract a frequency-smoothed version of the original spectrogram
# to use as filter
specEnv_orig = spectrogram(sound_orig, blur = c(300, 50),
 samplingRate = samplingRate, output = 'original', plot = TRUE)

# Synthesize source only, with flat spectrum
sound_unfilt = soundgen(sylLen = 2500, pitch = pitch,
  rolloff = 0, rolloffOct = 0,
  temperature = 0, formants = NULL, lipRad = 0,
  samplingRate = samplingRate,
  invalidArgAction = 'ignore')  # prevent soundgen from increasing samplingRate
# playme(sound_unfilt, samplingRate)
# meanSpectrum(sound_unfilt, samplingRate, yScale = 'max0')  # ~flat

# Force spectral envelope to the shape of target
sound_filt = addFormants(sound_unfilt, formants = NULL,
  formantFilter = specEnv_orig, samplingRate = samplingRate)
# playme(sound_filt, samplingRate)  # playme(sound_orig, samplingRate)
# spectrogram(sound_filt, samplingRate)  # spectrogram(sound_orig, samplingRate)

# The spectral envelope is now similar to the original recording. Compare:
par(mfrow = c(1, 2))
meanSpectrum(sound_orig, samplingRate, yScale = 'max0', alim = c(-50, 20))
meanSpectrum(sound_filt, samplingRate, yScale = 'max0', alim = c(-50, 20))
par(mfrow = c(1, 1))

## End(Not run)

Acoustic analysis

Description

Acoustic analysis of one or more sounds: pitch tracking, basic spectral characteristics, formants, estimated loudness (see getLoudness), roughness (see modulationSpectrum), novelty (see ssm), etc. The default values of arguments are optimized for human non-linguistic vocalizations. For high-precision work, first extract and manually correct pitch contours with pitch_app, PRAAT, or whatever, and then run analyze(pitchManual = ...) with these manual contours. For more information, see https://cogsci.se/soundgen/acoustic_analysis.html

Usage

analyze(
  x,
  samplingRate = NULL,
  scale = NULL,
  from = NULL,
  to = NULL,
  dynamicRange = 80,
  silence = 0.04,
  windowLength = 50,
  step = NULL,
  overlap = 50,
  wn = "gaussian",
  zp = 0,
  cutFreq = NULL,
  nFormants = 3,
  formants = list(),
  loudness = list(SPL_measured = 70),
  roughness = list(msType = "1D", specMethod = "spectrum", amRes = 1, specFun_pars =
    list(windowLength = 25, step = 2)),
  novelty = list(specFun = "melspec", kernelLen = 1000),
  pitchMethods = c("dom", "autocor"),
  pitchManual = NULL,
  pitchFloor = 75,
  pitchCeiling = 1000,
  priorMean = 300,
  priorSD = 6,
  priorAdapt = TRUE,
  nCands = 1,
  minVoicedCands = NULL,
  pitchDom = list(domThres = 0.1, domSmooth = 220),
  pitchAutocor = list(autocorThres = 0.7, autocorSmooth = 7, autocorUpsample = 25,
    autocorBestPeak = 0.975, interpol = "sinc"),
  pitchCep = list(cepThres = 0.75, cepZp = 0),
  pitchSpec = list(specThres = 0.05, specPeak = 0.25, specHNRslope = 0.8, specSmooth =
    150, specMerge = 0.1, specSinglePeakCert = 0.4, specRatios = 3),
  pitchHps = list(hpsNum = 5, hpsThres = 0.1, hpsNorm = 2, hpsPenalty = 2),
  pitchZc = list(zcThres = 0.1, zcWin = 5),
  harmHeight = list(harmThres = 3, harmTol = 0.25, harmPerSel = 5),
  subh = list(method = c("cep", "pitchCands", "harm")[1], nSubh = 5, tol = 0.05, nHarm =
    5, harmThres = 12, harmTol = 0.25),
  flux = list(thres = 0.15),
  amRange = c(10, 60),
  fmRange = NULL,
  shortestSyl = 20,
  shortestPause = 60,
  interpolPitch = list(win = 75, tol = 0.3, cert = 0.3),
  certWeight = 0.5,
  smooth = 1,
  smoothVars = c("pitch", "dom"),
  summaryFun = c("mean", "median", "sd"),
  invalidArgAction = c("adjust", "abort", "ignore"),
  reportEvery = NULL,
  cores = 1,
  plot = FALSE,
  osc = c("linear", "dB", "none"),
  showLegend = TRUE,
  savePlots = FALSE,
  embed = FALSE,
  pitchPlot = list(col = rgb(0, 0, 1, 0.75), lwd = 3, showPrior = TRUE),
  extraContour = NULL,
  ylim = NULL,
  xlab = "Time",
  ylab = NULL,
  main = NULL,
  width = 900,
  height = 500,
  units = "px",
  res = NA,
  ...
)

Arguments

Details

Each pitch tracker is controlled by its own list of settings, as follows:

pitchDom (lowest dominant frequency band)

• domThres (0 to 1) to find the lowest dominant frequency band, we do short-term FFT and take the lowest frequency with amplitude at least domThres

• domSmooth the width of smoothing interval (Hz) for finding dom

pitchAutocor (autocorrelation)

• autocorThres voicing threshold (unitless, ~0 to 1)

• autocorSmooth the width of smoothing interval (in bins) for finding peaks in the autocorrelation function. Defaults to 7 for sampling rate 44100 and smaller odd numbers for lower values of sampling rate

• autocorUpsample upsamples acf to this resolution (Hz) to improve accuracy in high frequencies

• autocorBestPeak amplitude of the lowest best candidate relative to the absolute max of the acf

• interpol method of interpolating the ACF: "sinc" for maximum precision, "none" for speed

pitchCep (cepstrum)

• cepThres voicing threshold (unitless, ~0 to 1)

• cepZp zero-padding of the spectrum used for cepstral pitch detection (final length of spectrum after zero-padding in points, e.g. 2 ^ 13)

pitchSpec (ratio of harmonics - BaNa algorithm)

• specThres voicing threshold (unitless, ~0 to 1)

• specPeak,specHNRslope when looking for putative harmonics in the spectrum, the threshold for peak detection is calculated as specPeak * (1 - HNR * specHNRslope)

• specSmooth the width of window for detecting peaks in the spectrum, Hz

• specMerge pitch candidates within specMerge semitones are merged with boosted certainty

• specSinglePeakCert (0 to 1) if F0 is calculated based on a single harmonic ratio (as opposed to several ratios converging on the same candidate), its certainty is taken to be specSinglePeakCert

• specMethod "commonFactor" = highest common factor of putative harmonics, "BaNa" = ratio of putative harmonics

• specRatios for method = "commonFactor", the number of harmonics and integer fractions to consider

pitchHps (harmonic product spectrum)

• hpsNum the number of times to downsample the spectrum

• hpsThres voicing threshold (unitless, ~0 to 1)

• hpsNorm the amount of inflation of hps pitch certainty (0 = none)

• hpsPenalty the amount of penalizing hps candidates in low frequencies (0 = none)

pitchZc (zero crossings)

• zcThres pitch candidates with certainty below this value are treated as noise and set to NA (0 = nothing discarded, 1 = pitch must be perfectly stable over zcWin)

• zcWin certainty in pitch candidates depends on how stable pitch is over zcWin glottal cycles (odd integer > 3)

Each of these lists also accepts graphical parameters that affect how pitch candidates are plotted, eg pitchDom = list(domThres = .5, col = 'yellow').

Other arguments that are lists of subroutine-specific settings include:

harmHeight (finding how high harmonics reach in the spectrum)

• harmThres minimum height of spectral peak, dB

• harmPerSel the number of harmonics per sliding selection

• harmTol maximum tolerated deviation of peak frequency from multiples of f0, proportion of f0

Value

A list with $detailed frame-by-frame descriptives and a $summary with one row per file, as determined by summaryFun (e.g., mean / median / SD of each acoustic variable across all STFT frames). Output measures include:

duration

total duration, s

duration_noSilence

duration from the beginning of the first non-silent STFT frame to the end of the last non-silent STFT frame, s (NB: depends strongly on windowLength and silence settings)

time

time of the middle of each frame (ms)

amEnvFreq,amEnvDep,amEnvPurity

frequency (Hz), purity (0 to 1), and depth (0 to 100) of amplitude modulation estimated from a smoothed amplitude envelope

amMsFreq,amMsPurity

frequency (Hz) and purity (dB) of amplitude modulation estimated via modulationSpectrum

ampl

root mean square of amplitude per frame, calculated as sqrt(mean(frame ^ 2))

ampl_noSilence

same as ampl, but ignoring silent frames

CPP

Cepstral Peak Prominence, dB (a measure of pitch quality, the ratio of the highest peak in the cepstrum to the regression line drawn through it)

dom

lowest dominant frequency band (Hz) (see "Pitch tracking methods / Dominant frequency" in the vignette)

entropyW

Wiener entropy of the spectrum of the current frame (=spectral flatness). Close to 0: pure tone or tonal sound with nearly all energy in harmonics; close to 1: white noise

entropySh

Normalized Shannon entropy of the spectrum of the current frame: 0 = pure tone, 1 = white noise

f1_freq, f1_width, ...

the frequency and bandwidth of the first nFormants formants per STFT frame, as calculated by phonTools::findformants

fluctuation

strength of low-frequency modulation at ~4 Hz (0.25-30 Hz), calculated from a modulation spectrum as a complement to psychoacoustic roughness; see modulationSpectrum

flux

feature-based flux, the rate of change in acoustic features such as pitch, HNR, etc. (0 = none, 1 = max); "epoch" is an audio segment between two peaks of flux that exceed a threshold of flux = list(thres = ...) (listed in output$detailed only)

fmFreq

frequency of frequency modulation (FM) such as vibrato or jitter, Hz

fmDep

depth of FM, semitones

fmPurity

purity or dominance of the main FM frequency (fmFreq), 0 to 1

harmEnergy

the amount of energy in upper harmonics, namely the ratio of total spectral mass above 1.25 x F0 to the total spectral mass below 1.25 x F0 (dB)

harmHeight

how high harmonics reach in the spectrum, based on the best guess at pitch (or the manually provided pitch values)

HNR

harmonics-to-noise ratio (dB), a measure of harmonicity (see "Pitch tracking methods / Autocorrelation"). If HNR = 0 dB, there is as much energy in harmonics as in noise

loudness

subjective loudness, in sone, corresponding to the chosen SPL_measured - see getLoudness

novelty

spectral novelty - a measure of how variable the spectrum is on a particular time scale, as estimated by ssm

peakFreq

the frequency with maximum spectral power (Hz)

pitch

post-processed pitch contour based on all F0 estimates

quartile25, quartile50, quartile75

the 25th, 50th, and 75th quantiles of the spectrum of voiced frames (Hz)

roughness

the amount of amplitude modulation in the roughness range, see modulationSpectrum and Anikin 2025

sharpness

psychoacoustic sharpness: related to spectral centroid, but calculated from a psychoacoustic loudness model, see getLoudness

specCentroid

the center of gravity of the frame's spectrum, first spectral moment (Hz)

specSlope

the slope of linear regression fit to the spectrum below cutFreq (dB/kHz)

subDep

estimated depth of subharmonics per frame: 0 = none, 1 = as strong as f0. NB: this depends critically on accurate pitch tracking

subRatio

the ratio of f0 to subharmonics frequency with strength subDep: 2 = period doubling, 3 = f0 / 3, etc.

voiced

is the current STFT frame voiced? TRUE / FALSE

References

Anikin, A. (2025) Acoustic estimation of voice roughness. Attention, Perception, & Psychophysics 87: 1771–1787.

See Also

pitch_app getLoudness segment getRMS

Examples

# Detailed documentation: https://cogsci.se/soundgen/acoustic_analysis.html

sound = soundgen(sylLen = 300, pitch = c(500, 400, 600),
  noise = list(time = c(0, 300), value = c(-40, 0)),
  temperature = 0.001,
  addSilence = 50)  # NB: always have some silence before and after!!!
# playme(sound, 16000)
a = analyze(sound, samplingRate = 16000, plot = TRUE)
str(a$detailed)  # frame-by-frame
a$summary        # summary per sound

## Not run: 
# For maximum processing speed (just basic spectral descriptives):
a = analyze(sound, samplingRate = 16000,
  plot = FALSE,         # no plotting
  pitchMethods = NULL,  # no pitch tracking
  loudness = NULL,      # no loudness analysis
  novelty = NULL,       # no novelty analysis
  roughness = NULL,     # no roughness analysis
  nFormants = 0         # no formant analysis
)

# Fancy plotting options:
a = analyze(sound, samplingRate = 44100, plot = TRUE,
  xlab = 'Time, ms', colorTheme = 'seewave', yScale = 'ERB',
  contrast = .5, ylim = c(0.05, 8), main = 'My plot',
  pitchMethods = c('dom', 'autocor', 'spec', 'hps', 'cep'),
  priorMean = NA,  # no prior info at all
  pitchDom = list(col = 'red', domThres = .25),
  pitchPlot = list(col = 'black', pch = 9, lty = 3, lwd = 3),
  extraContour = list(x = 'peakFreq', type = 'b', pch = 4, col = 'brown'),
  osc = 'dB', heights = c(2, 1))

# Analyze an entire folder in one go, saving spectrograms with pitch contours
# plus an html file for easy access
s2 = analyze('~/Downloads/temp',
  savePlots = TRUE,  # save the spectrograms with pitch contours
  showLegend = TRUE, yScale = 'bark',
  width = 20, height = 12,
  units = 'cm', res = 300, ylim = c(0, 5),
  cores = 4)  # use multiple cores to speed up processing
s2$summary[, 1:5]

# Analyzing ultrasounds (slow but possible, just adjust pitchCeiling)
s = soundgen(sylLen = 100, addSilence = 10,
  pitch = c(25000, 35000, 30000),
  formants = NA, rolloff = -12, rolloffKHz = 0,
  pitchSamplingRate = 350000, samplingRate = 350000, windowLength = 5,
  pitchCeiling = 45000, invalidArgAction = 'ignore',
  plot = TRUE)
# s is a bat-like ultrasound inaudible to humans

a = analyze(
  s, 350000, plot = TRUE,
  pitchFloor = 10000, pitchCeiling = 90000, priorMean = NA,
  pitchMethods = c('autocor', 'spec'),
  # probably shouldn't use pitchMethods = "dom" b/c of likely low-freq noise
  windowLength = 5, step = 2.5,
  shortestSyl = 10, shortestPause = 10,  # again, very short sounds
  interpolPitch = list(win = 10),  # again, very short sounds
  smooth = 0.1,  # might need less smoothing if very rapid f0 changes
  nFormants = 0, loudness = NULL, roughness = NULL, novelty = NULL)
# NB: ignore formants and loudness estimates for such non-human sounds

## End(Not run)

Annotation app

Description

Starts a shiny app for annotating audio. This is a simplified and faster version of formant_app intended only for making annotations. Use hotkeys (eg spacebar to play/stop) and avoid working with very large files.

Usage

annotation_app(...)

Arguments

Details

Recovering lost data

Every time you add an annotation or move in between files in the queue, the output you've got so far is saved in a temporary backup file. If the app crashes or is closed without saving the results, this backup file preserves your data. To recover it, restart annotation_app() - a dialog box will pop up and ask whether you want to append the old data to the new one. Even so, save your data regularly to be on the safe side!

Value

A list with two elements: $settings (a list of the last used settings) and $annotations (a data.frame with the annotations). Every time a new annotation is added, the app creates a backup csv file in the session's temporary directory, and it returns the final payload upon closing the app.

See Also

formant_app

Examples

## Not run: 
ann = annotation_app()  # runs in default browser such as Firefox or Chrome

ann = annotation_app(specType = 'reassigned', windowLength = 5, step = 1)

# full list of parameters that can be passed to annotation_app():
paste0(c(rownames(soundgen:::def_ann),
 'specType', 'spec_colorTheme', 'spec_yScale', 'osc', 'wn', 'audioMethod'),
 collapse = ', ')

# save the complete output, including the settings used
saveRDS(ann, 'my_annotations.rds')

# re-use the same settings in a future session
ann2 = do.call(annotation_app, ann$settings)

# To change system default browser, run something like:
options('browser' = '/usr/bin/firefox')  # path to the executable on Linux

## End(Not run)

Auditory spectrogram

Description

Produces an auditory spectrogram by convolving the sound with a bank of bandpass filters. The main difference from STFT is that we don't window the signal and de facto get variable temporal resolution in different frequency channels, as with a wavelet transform. The key settings are filterType, nFilters_oct, and yScale, which determine the type, number, and spacing of the filters, respectively. Gammatone filters were designed as a simple approximation of human perception - see Slaney 1993 "An Efficient Implementation of the Patterson–Holdsworth Auditory Filter Bank". Butterworth or Chebyshev filters are not meant to model perception, but can be useful for quickly plotting a sound.

Usage

audSpectrogram(
  x,
  samplingRate = NULL,
  scale = NULL,
  from = NULL,
  to = NULL,
  step = 10,
  dynamicRange = 80,
  filterType = c("gammatone", "butterworth", "chebyshev"),
  envelope = c("rms", "hil"),
  nFilters_oct = 6,
  nFilters = NULL,
  yScale = c("ERB", "bark", "mel", "log"),
  filterOrder = NULL,
  bandwidth = NULL,
  bandwidthMult = 1,
  minFreq = 20,
  maxFreq = NULL,
  minBandwidth = 10,
  output = c("all", "audSpec", "audSpec_processed", "filterbank", "filterbank_env",
    "filters"),
  reportEvery = NULL,
  cores = 1,
  plot = TRUE,
  savePlots = FALSE,
  embed = FALSE,
  plotFilters = FALSE,
  osc = c("linear", "dB", "none"),
  heights = c(3, 1),
  ylim = NULL,
  contrast = 0,
  brightness = 0,
  maxPoints = c(1e+05, 5e+05),
  colorTheme = "bw",
  col = NULL,
  extraContour = NULL,
  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

Value

A list for each analyzed file, including:

audSpec

auditory spectrogram: a matrix with frequency in rows (kHz) and time in columns (ms), offset by step/2

audSpec_processed

same dimensions, rescaled for plotting (log-transformed, contrast/brightness-adjusted, range 0–1)

filterbank

raw filter outputs: a matrix with one row per filter (ordered by center frequency) and one column per audio sample

filterbank_env

Hilbert envelopes of the filterbank, same dimensions as filterbank; NA if envelope = "rms"

filters

a dataframe giving the center frequencies, bandwidths, and lower/upper bounds of the used filters, all in Hz

Examples

data('speechEx', package = 'soundgen')

# auditory spectrogram
asp = audSpectrogram(speechEx, to = 1, step = 5)
dim(asp$audSpec)

# compare to STFT with similar time and frequency resolution (~100 times faster)
fs = spectrogram(speechEx, to = 1, yScale = 'ERB', windowLength = 5, step = 5)
dim(fs)

## Not run: 
# add bells and whistles
audSpectrogram(speechEx,
  nFilters = 128,
  dynamicRange = 150,
  osc = 'none',
  heights = c(2, 1),  # spectro/osc height ratio
  contrast = .4,  # increase contrast
  brightness = -.2,  # reduce brightness
  colorTheme = 'matlab',  # pick color theme...
  # col = hcl.colors(100, palette = 'Plasma'),  # ...or specify the colors
  cex.lab = .75, cex.axis = .75,  # text size and other base graphics pars
  grid = 5,  # to customize, add manually with graphics::grid()
  ylim = c(0.05, 8),  # always in kHz
  main = 'My auditory spectrogram' # title
  # + axis labels, etc
)

# NB: frequency resolution is controlled by both nFilters and bandwidth
audSpectrogram(speechEx, to = 1, nFilters = 15, bandwidth = 1/2)
audSpectrogram(speechEx, to = 1, nFilters = 15, bandwidth = 1/10)
audSpectrogram(speechEx, to = 1, nFilters = 100, bandwidth = 1/2)
audSpectrogram(speechEx, to = 1, nFilters = 100, bandwidth = 1/10)
audSpectrogram(speechEx, to = 1, nFilters_oct = 5, bandwidth = 1/10)
audSpectrogram(speechEx, to = 1, nFilters = 200, bandwidthMult = 1/3)

# caution: if bandwidths are too narrow relative to nFilters, there may be gaps
audSpectrogram(speechEx, to = 1, nFilters = 30, bandwidthMult = 1/3,
  plotFilters = TRUE, plot = FALSE)

# different filter types
audSpectrogram(speechEx, to = 1, filterType = 'gammatone')
audSpectrogram(speechEx, to = 1, filterType = 'butterworth')
audSpectrogram(speechEx, to = 1, filterType = 'chebyshev')

# save auditory spectrograms of all audio files in a folder
audSpectrogram('~/Downloads/temp', savePlots = TRUE, cores = 4)

## End(Not run)

Bandpass/stop filters

Description

Filtering in the frequency domain with FFT-iFFT: low-pass, high-pass, bandpass, and bandstop filters with flexible filter depth (dB) and adjustable bandwidths at filter boundaries. Similar to the "ffilter" function from the seewave package, but here we use FFT instead of STFT - that is, the entire sound is processed at once. This works best for relatively short sounds (seconds), but gives us maximum precision (e.g., for precise notch filtering) and doesn't affect the attack and decay. NAs are accepted and can be interpolated or preserved in the output. Because we don't do STFT, arbitrarily short vectors are also fine as input - for example, we can apply a low-pass filter prior to decimation when changing the sampling rate without aliasing. Note that, unlike pitchSmoothPraat, bandpass by default applies an abrupt cutoff instead of a smooth gaussian filter, but this behavior can be adjusted with the bw argument.

Usage

bandpass(
  x,
  samplingRate = NULL,
  lwr = NULL,
  upr = NULL,
  action = c("pass", "stop"),
  dB = Inf,
  bw = 0,
  na.rm = TRUE,
  from = NULL,
  to = NULL,
  normalize = FALSE,
  reportEvery = NULL,
  cores = 1,
  saveAudio = FALSE,
  plot = FALSE,
  savePlots = FALSE,
  embed = FALSE,
  width = 900,
  height = 500,
  units = "px",
  res = NA,
  ...
)

Arguments

Details

Algorithm: fill in NAs with constant interpolation at the edges and linear interpolation in the middle; perform FFT; set the frequency ranges to be filtered out to 0; perform inverse FFT; set to the original scale; put the NAs back in.

Value

Numeric vector of the same length as input.

Examples

# Filter white noise
s1 = fade(c(rnorm(2000, 0, 1), NA), samplingRate = 16000)

# low-pass
s1_lp = bandpass(s1, 16000, upr = 2000, bw = 100, plot = TRUE)
meanSpectrum(s1_lp, 16000)

# high-pass by 40 dB
s1_hp = bandpass(s1, 16000, lwr = 2000, dB = 40, bw = 50,
  plot = TRUE, windowLength = 5)
# windowLength is passed to meanSpectrum for plotting (ignore the warnings)

# bandstop
bandpass(s1, 16000, lwr = 1000, upr = 1800, bw = 50, action = 'stop', plot = TRUE)

# bandpass (turns noise into a whistle)
s1_bp = bandpass(s1, 16000, lwr = 2000, upr = 2100, bw = 25, plot = TRUE)
# playme(rep(s1_bp, 5))
# spectrogram(s1_bp, 16000)

# low-pass and interpolate a short vector with some NAs
x = rnorm(150, 10) + 3 * sin((1:50) / 5)
x[sample(seq_along(x), 50)] = NA
plot(x, type = 'l')
x_bandp = bandpass(x, samplingRate = 100, upr = 10, bw = 1)
points(x_bandp, type = 'l', col = 'blue')

# boost or attenuate a particular frequency band by 20 dB
s3 = bandpass(s1, 16000, lwr = 1700, upr = 2100, bw = 200,
  dB = 20, plot = TRUE)
# spectrogram(s3, 16000)
s4 = bandpass(s1, 16000, lwr = 2000, upr = 4300, bw = c(100, 500),
  dB = 60, action = 'stop', plot = TRUE)
# spectrogram(s4, 16000)

# precise notch filtering is possible, even in low frequencies
whiteNoise = rnorm(16000)
s5 = bandpass(whiteNoise, 16000, lwr = 30, upr = 40, normalize = TRUE,
              plot = TRUE, xlim = c(0, 0.5))
spectrum(s5, 16000, xlim = c(0, 0.1))
# playme(s5)
# spectrogram(s5, 16000, windowLength = 150, yScale = 'log')

## Not run: 
# high-pass all audio files in a folder
bandpass('~/Downloads/temp', lwr = 2000, saveAudio = TRUE, savePlots = TRUE)

## End(Not run)

Generate beat

Description

Generates percussive sounds from clicks through drum-like beats to sliding tones. The principle is to create a sine wave with rapid frequency modulation and to add a fade-out. No extra harmonics or formants are added. For this specific purpose, this is vastly faster and easier than to tinker with soundgen settings, especially since percussive syllables tend to be very short.

Usage

beat(
  nSyl = 10,
  sylLen = 200,
  pauseLen = 50,
  pitch = c(200, 10),
  samplingRate = 16000,
  fadeOut = TRUE,
  play = FALSE
)

Arguments

Value

The synthesized waveform as a numeric vector.

See Also

soundgen

Examples

playback = c(TRUE, FALSE)[2]
# a drum-like sound
s = beat(nSyl = 1, sylLen = 200,
         pitch = c(200, 100), play = playback)
# plot(s, type = 'l')

# a dry, muted drum
s = beat(nSyl = 1, sylLen = 200,
         pitch = c(200, 10), play = playback)

# sci-fi laser guns
s = beat(nSyl = 3, sylLen = 300,
         pitch = c(1000, 50), play = playback)

# machine guns
s = beat(nSyl = 10, sylLen = 10, pauseLen = 50,
         pitch = c(2300, 300), play = playback)

Clump a sequence into large segments

Description

clumper makes sure each homogeneous segment in a sequence is at least minLength long. Called by getIntegerRandomWalk(), addSubh(), naiveBayes(), etc. Algorithm: find the epochs shorter than minLength, merge max 1/4 of them with the largest neighbor, and repeat recursively until all epochs are at least minLength long. minLength can be a vector, in which case it is assumed to change over time.

Usage

clumper(x, minLength, n = length(x))

clumper2(x, minLength, n = length(x))

Arguments

Value

Returns the original sequence x transformed to homogeneous segments of required length, with the original class (e.g. character or factor).

Examples

s = c(1,3,2,2,2,0,0,4,4,1,1,1,1,1,3,3)
soundgen:::clumper(s, 2)
soundgen:::clumper(s, 3)
soundgen:::clumper(1:5, 10)
soundgen:::clumper(c('a','a','a','b','b','c','c','c','a','c'), 3)
soundgen:::clumper(x = c(1,2,1,2,1,1,1,1,3,1), minLength = c(1, 1, 1, 3))
soundgen:::clumper(as.factor(c('A','B','B','C')), 2)

## clumper vs clumper2
# clumper2 is faster for wiggly inputs
s = rpois(1000, 2)
plot(s, type = 'l')
plot(soundgen:::clumper(s, minLength = c(2, 3, 10)), type = 'l')
plot(soundgen:::clumper2(s, minLength = c(2, 3, 10)), type = 'l')
system.time(for (i in 1:10) soundgen:::clumper(s, minLength = c(2, 3, 10)))
system.time(for (i in 1:10) soundgen:::clumper2(s, minLength = c(2, 3, 10)))

# clumper1 is faster for step-like inputs
s2 = approx(c(2,4,3,2,0), n = 1001, method = 'constant')$y[1:1000]
plot(s2, type = 'l')
plot(soundgen:::clumper(s2, minLength = c(2, 3, 10)), type = 'l')
plot(soundgen:::clumper2(s2, minLength = c(2, 3, 10)), type = 'l')
system.time(for (i in 1:10) soundgen:::clumper(s2, minLength = c(2, 3, 10)))
system.time(for (i in 1:10) soundgen:::clumper2(s2, minLength = c(2, 3, 10)))

Compare sounds

Description

Computes distances between sounds based on comparing their spectrogram-like representations. compareSounds takes two sounds or feature matrices as input, whereas compareFolder takes a path to a folder with audio files or a list of feature matrices and returns a matrix of pairwise distances between them. Feature matrices are normalized and compared with Dynamic Time Warp (DTW), correlation, cosine distance, or pixel by pixel.

Usage

compareSounds(
  x,
  y,
  samplingRate = NULL,
  specFun = "melspec",
  specFun_pars = list(),
  logSpec = FALSE,
  method = c("cor", "cosine", "diff", "dtw"),
  padWith = NA,
  padDir = c("central", "left", "right"),
  dtw_pars = list()
)

compareFolder(
  myfolder = NULL,
  spectrograms = NULL,
  matchAllLengths = FALSE,
  specFun = "melspec",
  specFun_pars = list(),
  logSpec = FALSE,
  method = c("cor", "cosine", "diff", "dtw"),
  padWith = NA,
  padDir = c("central", "left", "right"),
  dtw_pars = list(),
  cores = 1,
  reportEvery = NULL
)

Arguments

Details

If the input is audio, several methods of producing spectrograms are available ("specFun"). For more customized options, just prepare your spectrograms or feature matrices first (time in columns, features like pitch, peak frequency, etc. in rows), and then pass them to compareSounds (see examples). All methods except for DTW require that the compared matrices should be of the same size. Compared sounds should ideally have the same sampling rate. If they differ, row (frequency bin) truncation is performed by position, keeping only the first min(nrow) rows, which approximates keeping frequencies up to the lower Nyquist frequency when both spectrograms use the same frequency resolution. In case of differences in duration, the shorter sound is padded with 0 (silence) or NA, as controlled by arguments padWith, padDir. If passing custom feature matrices, ensure that they have the same dimensions or that padding with 0 (silence) makes sense.

Value

compareSounds returns a dataframe with two columns: "method" for the method(s) used, and "distance" for the distance between the two sounds calculated with that method. The range of distances is [0, 1] for "cor", "cosine", and "diff", and [0, Inf) for "dtw". compareFolder returns a list of distance matrices (dist objects), one for each method.

Examples

s1 = soundgen(sylLen = 100, pitch = c(80, 180), formants = 'a')
s2 = soundgen(sylLen = 120, pitch = c(150, 350), formants = 'u')
compareSounds(s1, s2, samplingRate = 16000, method = c('cor', 'cosine', 'diff'))
# spectrogram(s1); playme(s1)
# spectrogram(s2); playme(s2)

## Not run: 
# NB: install the "dtw" library to run the examples

# compare all sounds in a folder, e.g.:
target = '~/Documents/Research/zz_test_audio/temp_long'
cf = compareFolder(target, logSpec = TRUE)
mds = as.data.frame(cmdscale(cf$cor))
plot(mds, type = 'n'); text(mds, labels = abbreviate(rownames(mds)))

# or use manually produced spectrograms
sp = spectrogram(target, windowLength = c(10, 40), overlap = 75,
  yScale = 'ERB', output = 'processed', plot = FALSE, cores = 4)
image(sp[[1]])
cf1 = compareFolder(spectrograms = sp)
mds1 = as.data.frame(cmdscale(cf1$cor))
plot(mds1, type = 'n'); text(mds1, labels = abbreviate(rownames(mds1)))

# extract a spectrogram-like representation using a custom function
# (e.g., full-resolution analytic envelopes instead of downsampled RMS)
compareSounds(s1, s2, samplingRate = 16000,
  specFun = function(x) matrix(hilbert_approx(x)$envelope, nrow = 1))

# some more examples
s1 = soundgen(formants = 'a', play = TRUE)
s2 = soundgen(formants = 'ae', play = TRUE)
s3 = soundgen(formants = 'eae', sylLen = 700, play = TRUE)
s4 = runif(8000, -1, 1)  # white noise
compareSounds(s1, s2, samplingRate = 16000)
compareSounds(s1, s4, samplingRate = 16000)

# the central section of s3 is more similar to s1 than is the beg/end of s3
compareSounds(s1, s3, samplingRate = 16000, padDir = 'left')
compareSounds(s1, s3, samplingRate = 16000, padDir = 'central')

# padding with 0 penalizes differences in duration, whereas padding with NA
# is like saying we only care about the overlapping part
compareSounds(s1, s4, samplingRate = 16000, padWith = 0)
compareSounds(s1, s4, samplingRate = 16000, padWith = NA)

# different types of spectrograms produce quite different results
compareSounds(s1, s3, samplingRate = 16000, specFun = 'stft')
compareSounds(s1, s3, samplingRate = 16000, specFun = 'melspec')
compareSounds(s1, s3, samplingRate = 16000, specFun = 'mfcc')
compareSounds(s1, s3, samplingRate = 16000, specFun = 'audSpec')

# pass additional control parameters to specFun and DTW
compareSounds(s1, s3, samplingRate = 16000,
              specFun = 'melspec',
              specFun_pars = list(nbands = 128),
              dtw_pars = list(dist.method = "Manhattan"))

# use feature matrices instead of spectrograms
# (time in columns, features in rows)
a1 = t(as.matrix(analyze(s1, samplingRate = 16000)$detailed))
a1 = a1[4:nrow(a1), ]; a1[is.na(a1)] = 0  # don't use dur and time stamps
a2 = t(as.matrix(analyze(s2, samplingRate = 16000)$detailed))
a2 = a2[4:nrow(a2), ]; a2[is.na(a2)] = 0
a4 = t(as.matrix(analyze(s4, samplingRate = 16000)$detailed))
a4 = a4[4:nrow(a4), ]; a4[is.na(a4)] = 0
compareSounds(a1, a2, method = c('cosine', 'dtw'))
compareSounds(a1, a4, method = c('cosine', 'dtw'))

## End(Not run)

Join two waveforms by cross-fading

Description

crossFade joins two input vectors (waveforms) by overlapping and cross-fading them. This is standard for joining two separately generated sounds without audible artifacts.

Usage

crossFade(
  x,
  y,
  crossLen = 15,
  crossLen_points = NULL,
  samplingRate = NULL,
  shape = c("equalPower", "lin", "exp", "log", "cos", "logistic", "gaussian"),
  cutToZC = FALSE,
  steepness = 1,
  SNR = NULL
)

Arguments

Value

A numeric vector with the original sampling rate. The length is going to be length(x) + length(y) - crossLen_points, or a bit less if cutting to nearest zero crossings.

See Also

fade

Examples

sound1 = sin(1:100 / 9)
sound2 = sin(7:106 / 3)
plot(c(sound1, sound2), type = 'b')
# an ugly discontinuity at 100 that will make an audible click

sound = crossFade(sound1, sound2, crossLen_points = 25)
plot(sound, type = 'b') # a nice, smooth transition
length(sound) # but note that cross-fading costs us crossLen_points = 25
# points even without cutting to zc, and much more if we trim to
# zero crossings and then overlap

sound = crossFade(sound1, sound2, crossLen_points = 25,
  SNR = 6)  # make sound1 6 dB louder (RMS) than sound2
plot(sound, type = 'b') # a nice, smooth transition

# Actual sounds, alternative shapes of fade-in/out
sound3 = soundgen(formants = 'a', pitch = 200, sylLen = 300,
                  addSilence = 0, attackLen = c(50, 0))
sound4 = soundgen(formants = 'u', pitch = 200, sylLen = 300,
                  addSilence = 0, attackLen = c(0, 50))

# simple concatenation (with a click)
cf1 = c(sound3, sound4)
osc(cf1, 16000, from = .25, to = .35)
# playme(cf1, 16000)

# concatenation from zc to zc (no click, but a rough transition)
cf2 = crossFade(sound3, sound4, crossLen_points = 0, cutToZC = TRUE)
osc(cf2, 16000, from = .25, to = .35)
# playme(cf2, 16000)

# crossFade over 35 ms - brief, but smooth
cf3 = crossFade(sound3, sound4, crossLen = 35, samplingRate = 16000)
osc(cf3, 16000, from = .25, to = .35)
# playme(cf3, 16000)

# equal-power cross-fade over 300 ms (shortens the sound by ~300 ms)
cf4 = crossFade(sound3, sound4, samplingRate = 16000,
                 crossLen = 300, shape = 'equalPower')
osc(cf4, 16000)
# playme(cf4, 16000)

Detect NLP

Description

(Experimental) A function for automatically detecting and annotating nonlinear vocal phenomena (NLP). Algorithm: analyze the audio using analyze and phasegram, then use the extracted frame-by-frame descriptives to classify each frame as having no NLP ("none"), subharmonics ("sh"), sidebands / amplitude modulation ("sb"), or deterministic chaos ("chaos"). The classification is performed by a naiveBayes algorithm adapted to autocorrelated time series and pretrained on a manually annotated corpus of vocalizations. Whenever possible, check and correct pitch tracks prior to running the algorithm. See naiveBayes for tips on using adaptive priors and clumping to account for the fact that NLP typically occur in continuous segments spanning multiple frames.

Usage

detectNLP(
  x,
  samplingRate = NULL,
  predictors = NULL,
  thresProb = 0.4,
  voicelessToNone = FALSE,
  train = soundgen::detectNLP_training_nonv,
  scale = NULL,
  from = NULL,
  to = NULL,
  pitchManual = NULL,
  analyze_pars = list(windowLength = 50, roughness = list(msType = "1D", specMethod =
    "spectrum", specFun_pars = list(windowLength = 25, step = 2), amRes = NULL,
    roughRange = NULL, roughMean = 100, roughSD = 8), amRange = c(20, 150)),
  phasegram_pars = list(nonlinStats = "d2"),
  naiveBayes_pars = list(prior = "static", wlClumper = 3),
  jumpThres = 14,
  jumpWindow = 100,
  reportEvery = NULL,
  cores = 1,
  plot = FALSE,
  savePlots = FALSE,
  embed = FALSE,
  main = NULL,
  xlab = NULL,
  ylab = NULL,
  ylim = NULL,
  width = 900,
  height = 500,
  units = "px",
  res = NA,
  ...
)

Arguments

Value

A list of datasets, one per input file, with acoustic descriptives per frame (returned by analyze and phasegram), probabilities of each NLP type per frame, and the putative classification of NLP per frame. The time step is equal to the larger of the steps passed to analyze() and phasegram().

Examples

## Not run: 
target = soundgen(sylLen = 2000, addSilence = 0, temperature = 1e-2,
  pitch = c(380, 550, 500, 220), subDep = c(0, 0, 40, 0, 0, 0, 0, 0),
  amDep = c(0, 0, 0, 0, 80, 0, 0, 0), amFreq = 80,
  noise = c(-10, rep(-40, 5)),
  jitterDep = c(0, 0, 0, 0, 0, 3),
  plot = TRUE, play = TRUE)

# classifier trained on manually annotated recordings of human nonverbal
# vocalizations
nlp = detectNLP(target, 16000,
  plot = TRUE, ylim = c(0, 4))

# classifier trained on synthetic, soundgen()-generated sounds
nlp = detectNLP(target, 16000,
  train = soundgen::detectNLP_training_synth,
  plot = TRUE, ylim = c(0, 4))
head(nlp[, c('time', 'pr')])
table(nlp$pr)
plot(nlp$amEnvDep, type = 'l')
plot(nlp$none, type = 'l')
points(nlp$sb, type = 'l', col = 'blue')
points(nlp$sh, type = 'l', col = 'green')
points(nlp$chaos, type = 'l', col = 'red')

# detection of pitch jumps
s1 = soundgen(sylLen = 1200, temperature = .001, pitch = list(
  time = c(0, 350, 351, 890, 891, 1200),
  value = c(140, 230, 460, 330, 220, 200)))
playme(s1, 16000)
nlp1 = detectNLP(s1, 16000, plot = TRUE, ylim = c(0, 3))

## End(Not run)

Training datasets for NLP detection

Description

The results of running naiveBayes_train on acoustically analyzed 969 human nonverbal vocalizations (detectNLP_training_nonv) and on 5000 synthetic sounds with or without NLP created with soundgen (detectNLP_training_synth). Used by detectNLP.

Usage

detectNLP_training_nonv

detectNLP_training_synth

Format

An object of class list of length 12.

An object of class list of length 12.

Estimate vocal tract length

Description

Estimates the length of vocal tract based on formant frequencies. If method = 'meanFormant', vocal tract length (VTL) is calculated separately for each formant, and then the resulting VTLs are averaged. The equation used is (2 * n - 1) * speedSound / (4 * Fn) for a closed-open tube (mouth open) and n * speedSound / (2 * Fn) for an open-open or closed-closed tube (eg closed mouth in mmm or open mouth and open glottis in whispering). If method = 'meanDispersion', formant dispersion is calculated as the mean distance between formants, and then VTL is calculated as speed of sound / 2 / formant dispersion. If method = 'regression', formant dispersion is estimated using the regression method described in Reby et al. (2005) and Anikin et al. (2024). For a review of VTL-related summary measures of formant frequencies, refer to Pisanski et al. (2014). See also schwa for VTL estimation with additional information on formant frequencies.

Usage

estimateVTL(
  formants,
  method = c("regression", "meanDispersion", "meanFormant"),
  interceptZero = TRUE,
  tube = c("closed-open", "open-open", "open-closed", "closed-closed"),
  speedSound = 35400,
  checkFormat = TRUE,
  output = c("simple", "detailed"),
  plot = FALSE
)

Arguments

Value

If output = 'simple' (default), returns the estimated vocal tract length in cm. If output = 'detailed' and method = 'regression', returns a list with extra stats used for plotting. Namely, $regressionInfo$infl gives the influence of each observation calculated as the absolute change in VTL with vs without the observation * 10 + 1 (the size of labels on the first plot). $vtlPerFormant$vtl gives the VTL as it would be estimated if only the first nFormants were used.

References

• Reby, D., McComb, K., Cargnelutti, B., Darwin, C., Fitch, W. T., & Clutton-Brock, T. (2005). Red deer stags use formants as assessment cues during intrasexual agonistic interactions. Proceedings of the Royal Society B: Biological Sciences, 272(1566), 941-947.

• Pisanski et al. (2014). Vocal indicators of body size in men and women: a meta-analysis. Animal Behaviour, 95, 89-99.

• Anikin, A., Barreda, S. & Reby, D. (2024) A practical guide to calculating vocal tract length and scale-invariant formant patterns. Behavior Research Methods 56, 5588–5604.

See Also

schwa

Examples

estimateVTL(NA)
estimateVTL(500)
estimateVTL(c(600, 1850, 2800, 3600, 5000), plot = TRUE)
es = estimateVTL(c(600, 1850, 2800, 3600, 5000), plot = TRUE, output = 'detailed')
es

estimateVTL(c(1200, 2000, 2800, 3800, 5400, 6400),
  tube = 'open-open', interceptZero = FALSE, plot = TRUE)
estimateVTL(c(1200, 2000, 2800, 3800, 5400, 6400),
  tube = 'open-open', interceptZero = TRUE, plot = TRUE)

# Multiple measurements are OK
estimateVTL(
  formants = list(f1 = c(540, 600, 550),
  f2 = 1650, f3 = c(2400, 2550)),
  plot = TRUE, output = 'detailed')
# NB: this is better than averaging formant values. Cf.:
estimateVTL(
  formants = list(f1 = mean(c(540, 600, 550)),
  f2 = 1650, f3 = mean(c(2400, 2550))),
  plot = TRUE)

# Missing values are OK
estimateVTL(c(600, 1850, 3100, NA, 5000), plot = TRUE)
estimateVTL(list(f1 = 500, f2 = c(1650, NA, 1400), f3 = 2700), plot = TRUE)

# Note that VTL estimates based on the commonly reported 'meanDispersion'
# depend only on the first and last formants
estimateVTL(c(500, 1400, 2800, 4100), method = 'meanDispersion')
estimateVTL(c(500, 1100, 2300, 4100), method = 'meanDispersion') # identical
# ...but this is not the case for 'meanFormant' and 'regression' methods
estimateVTL(c(500, 1400, 2800, 4100), method = 'meanFormant')
estimateVTL(c(500, 1100, 2300, 4100), method = 'meanFormant') # much longer

# Compare the results produced by the three methods
nIter = 100
out = data.frame(meanFormant = rep(NA, nIter), meanDispersion = NA, regression = NA)
for (i in 1:nIter) {
  # generate a random formant configuration
  f = runif(1, 300, 900) + (1:6) * rnorm(6, 1000, 200)
  out$meanFormant[i]    = estimateVTL(f, method = 'meanFormant')
  out$meanDispersion[i] = estimateVTL(f, method = 'meanDispersion')
  out$regression[i]     = estimateVTL(f, method = 'regression')
}
pairs(out)
cor(out)
# 'meanDispersion' is pretty different, while 'meanFormant' and 'regression'
# give broadly comparable results

Fade

Description

Applies fade-in and/or fade-out of variable length, shape, and steepness. The resulting effect softens the attack and release of a waveform.

Usage

fade(
  x,
  fadeIn = 50,
  fadeOut = 50,
  samplingRate = NULL,
  scale = NULL,
  shape = c("lin", "exp", "log", "cos", "logistic", "gaussian"),
  steepness = 1,
  reportEvery = NULL,
  cores = 1,
  saveAudio = FALSE,
  plot = FALSE,
  savePlots = FALSE,
  embed = FALSE,
  width = 900,
  height = 500,
  units = "px",
  res = NA
)

Arguments

Value

A numeric vector of the same length as input

See Also

crossFade

Examples

# Fading a real sound: say we want fast attack and slow release
s = soundgen(attackLen = 0, windowLength = 10,
             sylLen = 500, addSilence = 0)
# playme(s)
s1 = fade(s, fadeIn = 40, fadeOut = 350,
          samplingRate = 16000, shape = 'cos', plot = TRUE)
# playme(s1)

# Illustration of fade shapes
x = runif(4000, min = -1, max = 1)  # make sure to zero-center input!!!
# plot(x, type = 'l')
y = fade(x, samplingRate = 1000, fadeIn = 1000, fadeOut = 0, plot = TRUE)
y = fade(x, samplingRate = 1000, fadeIn = 1000, fadeOut = 1500,
         shape = 'exp', steepness = 1, plot = TRUE)
y = fade(x, samplingRate = 1000, fadeIn = 1500, fadeOut = 500,
         shape = 'log', steepness = 1, plot = TRUE)
y = fade(x, samplingRate = 1000, fadeIn = 1500, fadeOut = 500,
         shape = 'log', steepness = 3, plot = TRUE)
y = fade(x, samplingRate = 1000, fadeIn = 1500, fadeOut = 1500,
         shape = 'cos', plot = TRUE)
y = fade(x, samplingRate = 1000, fadeIn = 1500, fadeOut = 1500,
         shape = 'logistic', steepness = 1, plot = TRUE)
y = fade(x, samplingRate = 1000, fadeIn = 1500, fadeOut = 1500,
         shape = 'logistic', steepness = 3, plot = TRUE)
y = fade(x, samplingRate = 1000, fadeIn = 1500, fadeOut = 1500,
         shape = 'gaussian', steepness = 1.5, plot = TRUE)

## Not run: 
fade('~/Downloads/temp', fadeIn = 500, fadeOut = 500, savePlots = TRUE)

## End(Not run)

Filter sound by modulation spectrum

Description

Manipulates the modulation spectrum (MS) of a sound to attenuate or amplify certain frequencies of amplitude modulation (AM) and frequency modulation (FM). Algorithm: produces a modulation spectrum with modulationSpectrum, modifies it with filterMS, converts the modified MS back to a spectrogram, and finally inverts the spectrogram with invertSpectrogram, thus producing a sound with (approximately) the desired characteristics of the MS. Note that the last step of inverting the spectrogram introduces some noise, so the resulting MS is not precisely the same as the intermediate filtered version. In practice this means that some residual energy will still be present in the filtered-out frequency range (see examples).

Usage

filterSoundByMS(
  x,
  samplingRate = NULL,
  from = NULL,
  to = NULL,
  logSpec = TRUE,
  windowLength = 25,
  step = NULL,
  overlap = 80,
  wn = "hanning",
  amCond = NULL,
  fmCond = NULL,
  jointCond = NULL,
  action = c("stop", "pass"),
  dynamicRange = 80,
  dB = Inf,
  smooth = NULL,
  protect = c("level", "meanSpectrum"),
  warnAsymmetric = TRUE,
  initialPhase = c("orig", "spsi", "random", "zero"),
  nIter = 50,
  reportEvery = NULL,
  cores = 1,
  play = FALSE,
  saveAudio = FALSE,
  plot = TRUE,
  savePlots = FALSE,
  embed = FALSE,
  width = 900,
  height = 500,
  units = "px",
  res = NA
)

filterMS(
  ms,
  amCond = NULL,
  fmCond = NULL,
  jointCond = NULL,
  action = c("stop", "pass"),
  dB = Inf,
  smooth = NULL,
  protect = c("level", "meanSpectrum"),
  warnAsymmetric = TRUE,
  plot = TRUE
)

Arguments

Details

By default, the filter mask is smoothed to avoid abrupt cutoffs. Smoothing is specified in physical units: FM in cycles/kHz and AM in Hz. Set smooth = NA to disable smoothing. If logSpec = TRUE, filtering is performed on the modulation spectrum of the log-magnitude spectrogram, which is often more robust and more physiologically plausible than filtering linear magnitudes. If initialPhase = 'orig', the phase of the original sound is used. With nIter = 0, the original phase is used directly; with nIter > 0, the original phase is used to initialize the Griffin-Lim iterations.

Value

Filtered audio as a numeric vector normalized to [-1, 1] with the same sampling rate as input.

Choosing filtering settings

MS filtering operates globally on the spectrogram's envelope, so any transient whose modulation content overlaps the filtered band — onsets/offsets, plosives, amplitude steps — is partially removed together with the target, producing a ripple that decays inward from the transient. Smoothing (blurring the filter mask before applying it) shortens the ripple's tails but cannot eliminate it completely. Two parameters control the trade-off between successful filtering and introducing spectro-temporal artifacts: (1) dynamicRange: lower values shrink the silence-to-sound step in the log domain and lift the linear target above the zero-clipping knee, suppressing edge ripple at the cost of a uniform residual noise floor, most visible in silent regions, and (2) dB: finite attenuation, e.g. 20–40 dB, removes proportionally less of the transients' in-band content than dB = Inf. Use dB = Inf with a large dynamicRange for steady, sustained sounds and finite dB with dynamicRange = 60–80 for sounds with transients, such as speech. Use logSpec = TRUE for manipulating the temporal envelope / AM in a perceptually meaningful manner. logSpec = FALSE weights components by absolute magnitude and avoids amplifying quiet regions, but shows stronger edge artifacts. Choose protect = c('level', 'meanSpectrum') (the default setting) when filtering AM. Do not protect 'envelope', which contains the very AM being removed. Set protect = c('level', 'envelope') when filtering FM, so as to keep the natural amplitude contour. Drop 'meanSpectrum' only if you also want to flatten the static spectral shape. Whenever the condition includes 0 on an axis (e.g. removing slow AM with abs(am) < 3), keep the corresponding protection to avoid a collapsed, "AC-coupled" reconstruction. For mild filtering, initialPhase = 'orig' with nIter = 0 is fast and preserves the original phase; increase nIter for stronger filtering.

See Also

invertSpectrogram

Examples

# Create a sound to be filtered
s = soundgen(sylLen = 500, samplingRate = 16000,
    amFreq = 25, amDep = 50,
    addSilence = 50, plot = TRUE)
# playme(s, 16000)

# Filter to remove the rapid AM at 25 Hz
s_filt = fade(filterSoundByMS(s, samplingRate = 16000,
  amCond = 'abs(am) > 22 & abs(am) < 28',
  action = 'stop',
  plot = TRUE))
# NB: clicks are often introduced at the onset and offset - fade and renormalize
s_filt = s_filt / max(abs(s_filt))
# playme(s_filt, samplingRate = 16000)
spectrogram(s_filt, 16000, windowLength = 25)

# Using filterMS()
ms = modulationSpectrum(s, 16000, returnComplex = TRUE)$detailed$complex

# Remove all AM over 25 Hz
filterMS(ms, amCond = 'abs(am) > 25', protect = FALSE)
filterMS(ms, amCond = 'abs(am) > 25', protect = TRUE)
filterMS(ms, amCond = 'abs(am) > 25', protect = TRUE, smooth = NA)

# amCond and fmCond are OR-conditions
filterMS(ms, amCond = 'abs(am) > 15', fmCond = 'abs(fm) > 5', action = 'stop')
filterMS(ms, amCond = 'abs(am) > 15', fmCond = 'abs(fm) > 5', action = 'pass')
filterMS(ms, amCond = 'abs(am) > 10 & abs(am) < 25', action = 'stop')

# jointCond is more flexible
filterMS(ms, jointCond = 'am * fm < 5', action = 'stop')
filterMS(ms, jointCond = 'am^2 + (fm*3)^2 < 200', action = 'pass')

# So:
filterMS(ms, jointCond = 'abs(am) > 5 | abs(fm) < 5')  # general
# ...is the same as:
filterMS(ms, amCond = 'abs(am) > 5', fmCond = 'abs(fm) < 5')  # faster

## Not run: 
data('speechEx', package = 'soundgen')
samplingRate = speechEx@samp.rate
playme(speechEx)
spectrogram(speechEx)

# Remove AM above 3 Hz from a bit of speech (removes most temporal details)
s_filt1 = fade(filterSoundByMS(speechEx, amCond = 'abs(am) > 3',
                          action = 'stop',
                          smooth = c(0.5, 2)))
spectrogram(s_filt1, samplingRate = samplingRate)
playme(s_filt1, samplingRate)

# Remove FM above 3 cycles/kHz (masks f0 but preserves formants)
s_filt3 = fade(filterSoundByMS(speechEx,
  fmCond = 'abs(fm) > 3', action = 'stop',
  logSpec = TRUE, smooth = c(1, 5),
  protect = c('level', 'envelope')))
spectrogram(s_filt3, samplingRate = samplingRate)
playme(s_filt3, samplingRate)

# Joint spectro-temporal filtering with a symmetric condition
s_filt4 = fade(filterSoundByMS(speechEx,
  jointCond = 'am^2 + (fm * 3)^2 < 300', action = 'stop',
  logSpec = TRUE, dynamicRange = 60, smooth = c(1, 10),
  protect = c('level', 'meanSpectrum')))
spectrogram(s_filt4, samplingRate = samplingRate)
playme(s_filt4, samplingRate)

# Process all files in a folder, save filtered audio and plots
s_filt = filterSoundByMS('~/Downloads/temp2',
  saveAudio = TRUE, savePlots = TRUE,
  amCond = 'abs(am) > 15', fmCond = 'abs(fm) > 5',
  action = 'stop', smooth = c(0.5, 5),
  nIter = 10)

## End(Not run)

Find peaks

Description

A bare-bones, very fast function to find local maxima (peaks) in a numeric vector. NAs are set to the minimum value of input and effectively ignored. The first and last values can also be peaks. Plateaus are NOT treated as peaks.

Usage

findPeaks(
  x,
  wl = 3,
  thres = NULL,
  thresRel = NULL,
  dropFirstLast = FALSE,
  plot = FALSE
)

Arguments

Value

A vector of indices of local maxima.

Examples

x = c(5, 2, 3, 3, 2, 7, 3, 2, 5, 5, 1, 0, 2, NA)
findPeaks(x, wl = 3, plot = TRUE)
findPeaks(x, wl = 3, thres = 2, plot = TRUE)
findPeaks(x, wl = 3, thres = 100, plot = TRUE)
findPeaks(x, wl = 4, plot = TRUE)
findPeaks(x, wl = 6, plot = TRUE)
findPeaks(x, wl = 3, thresRel = 2, plot = TRUE)
findPeaks(x, wl = 3, thresRel = 4, plot = TRUE)
findPeaks(x, wl = 3, thres = 6, thresRel = 2, plot = TRUE)
findPeaks(x, dropFirstLast = TRUE)

Flat envelope / compressor

Description

Applies a compressor - that is, flattens the amplitude envelope of a waveform, reducing the difference in amplitude between loud and quiet sections. This is achieved by dividing the waveform by some function of its smoothed amplitude envelope (Hilbert, peak or root mean square).

Usage

flatEnv(
  x,
  samplingRate = NULL,
  scale = NULL,
  compression = 1,
  method = c("hil", "rms", "peak"),
  windowLength = 50,
  killDC = FALSE,
  dynamicRange = 40,
  reportEvery = NULL,
  cores = 1,
  saveAudio = FALSE,
  plot = FALSE,
  savePlots = FALSE,
  embed = FALSE,
  col = "blue",
  width = 900,
  height = 500,
  units = "px",
  res = NA,
  ...
)

compressor(
  x,
  samplingRate = NULL,
  scale = NULL,
  compression = 1,
  method = c("hil", "rms", "peak"),
  windowLength = 50,
  killDC = FALSE,
  dynamicRange = 40,
  reportEvery = NULL,
  cores = 1,
  saveAudio = FALSE,
  plot = FALSE,
  savePlots = FALSE,
  embed = FALSE,
  col = "blue",
  width = 900,
  height = 500,
  units = "px",
  res = NA,
  ...
)

Arguments

Value

If the input is a single audio (file, Wave, or numeric vector), returns the compressed waveform as a numeric vector with the original sampling rate and scale. If the input is a folder with several audio files, returns a list of compressed waveforms, one for each file.

See Also

transplantEnv

Examples

a = rnorm(500) * seq(1, 0, length.out = 500)
b = flatEnv(a, 1000, plot = TRUE, windowLength = 5)    # too short
c = flatEnv(a, 1000, plot = TRUE, windowLength = 450)  # too long
d = flatEnv(a, 1000, plot = TRUE, windowLength = 100)  # about right

## Not run: 
s = soundgen(sylLen = 1000, ampl = c(0, -40, 0), plot = TRUE)
# playme(s)
s_flat1 = flatEnv(s, 16000, dynamicRange = 60, plot = TRUE,
                  windowLength = 50, method = 'hil')
s_flat2 = flatEnv(s, 16000, dynamicRange = 60, plot = TRUE,
                  windowLength = 10, method = 'rms')
s_flat3 = flatEnv(s, 16000, dynamicRange = 60, plot = TRUE,
                  windowLength = 10, method = 'peak')
# playme(s_flat2)

# Remove DC offset
s1 = c(rep(0, 50), runif(1000, -1, 1), rep(0, 50)) +
     seq(.3, 1, length.out = 1100)
s2 = flatEnv(s1, 16000, plot = TRUE, windowLength = 50, killDC = FALSE)
s3 = flatEnv(s1, 16000, plot = TRUE, windowLength = 50, killDC = TRUE)

# Compress and save all audio files in a folder
s4 = flatEnv('~/Downloads/temp',
             method = 'peak', compression = .5,
             saveAudio = TRUE,
             savePlots = TRUE,
             col = 'green', lwd = 5)
osc(s4[[1]])

## End(Not run)

Interactive formant tracker

Description

Starts a shiny app for manually correcting formant measurements. For more tips, see pitch_app and http://cogsci.se/soundgen.html.

Usage

formant_app(...)

Arguments

Details

Suggested workflow: load one or several audio files (wav/mp3), preferably not longer than a minute or so. Select a region of interest in the spectrogram - for example, a sustained vowel with clear and relatively steady formants. Double-click within the selection to create a new annotation (you may add a text label if needed). If you are satisfied with the automatically calculated formant frequencies, proceed to the next region of interest. If not, there are three ways to adjust them: (1) click the spectrogram within selection (pick the formant number to adjust by clicking the formant boxes); (2) single-click the spectrum to use the cursor's position, or (3) double-click the spectrum to use the nearest spectral peak. When done with a file, move on to the next one in the queue. Use the orange button to download the results. To continue work, upload the output file from the previous session together with the audio files (you can rename it, but keep the .csv extension). Use hotkeys (eg spacebar to play/stop) and avoid working with very large files.

Recovering lost data

Every time you add an annotation or move in between files in the queue, the output you've got so far is saved in a temporary backup file. If the app crashes or is closed without saving the results, this backup file preserves your data. To recover it, restart formant_app() - a dialog box will pop up and ask whether you want to append the old data to the new one. Even so, save your data regularly to be on the safe side!

Value

A list of the last used settings ($settings) plus a data.frame with the formant measurements. Every time a new annotation is added, the app creates a backup csv file in the session's temporary directory, and it returns the final payload upon closing the app.

See Also

pitch_app annotation_app

Examples

## Not run: 
f = formant_app()  # runs in default browser such as Firefox or Chrome

f1 = formant_app(specType = 'reassigned', windowLength = 5, step = 1)

# run the app with previously used settings
f2 = do.call(formant_app, f1$settings)

# save the complete output, including the settings used
saveRDS(f2, 'my_formant_analysis.rds')

# To change system default browser, run something like:
options('browser' = '/usr/bin/firefox')  # path to the executable on Linux

## End(Not run)

Gaussian smoothing in 2D

Description

Takes a matrix of numeric values and smooths it by convolution with a symmetric Gaussian window function. Values outside the matrix are either treated as zero, attenuating the edges, or assumed to continue beyond the edges.

Usage

gaussianSmooth2D(
  m,
  kernelSize = 5,
  kernelSD = 0.5,
  action = c("blur", "unblur"),
  amount = 0.5,
  padWith = c("repeat", "zero"),
  plotKernel = FALSE
)

Arguments

Value

A numeric matrix of the same dimensions as input.

See Also

modulationSpectrum spectrogram

Examples

data('speechEx', package = 'soundgen')
s = spectrogram(speechEx, from = 0, to = 1, windowLength = 10,
  output = 'original', plot = FALSE)
s = log(s + .001)
image(t(s))
s1 = gaussianSmooth2D(s, kernelSize = 5, plotKernel = TRUE)
image(t(s1))

# more smoothing in time than in frequency
s2 = gaussianSmooth2D(s, kernelSize = c(5, 15))
image(t(s2))

# vice versa - more smoothing in frequency
s3 = gaussianSmooth2D(s, kernelSize = c(25, 3))
image(t(s3))

# smoothing only in one dimension
s4 = gaussianSmooth2D(s, kernelSize = c(25, 0))
image(t(s4))
s5 = gaussianSmooth2D(s, kernelSize = c(0, 15))
image(t(s5))

# sharpen the image
s6 = gaussianSmooth2D(s, kernelSize = 5, action = 'unblur', amount = .5)
image(t(s6))

Generate noise

Description

Generates noise of length len and with spectrum defined by rolloff parameters OR by a specified filter formantFilter. This function is called internally by soundgen, but it may be more convenient to call it directly when synthesizing non-biological noises defined by specific spectral and amplitude envelopes rather than formants: the wind, whistles, impact noises, etc. See beat for similarly simplified functions for tonal non-biological sounds.

Usage

generateNoise(
  len,
  rolloffNoise = -4,
  noiseFlatSpec = 1200,
  rolloffNoiseExp = 0,
  formantFilter = NULL,
  noise = NULL,
  attackLen = 10,
  samplingRate = 16000,
  windowLength = 50,
  step = NULL,
  overlap = 75,
  wn = "gaussian",
  smoothing = list(),
  play = FALSE
)

Arguments

Details

Algorithm: paints a spectrogram with desired characteristics, sets phase to zero, and generates a time sequence via inverse FFT.

Value

The generated waveform as a numeric vector.

See Also

soundgen beat

Examples

# .5 s of white noise
samplingRate = 16000
noise1 = soundgen:::generateNoise(len = samplingRate * .5,
  samplingRate = samplingRate)
meanSpectrum(noise1, samplingRate)
# playme(noise1, samplingRate)

# Percussion (run a few times to notice stochasticity due to temperature = .25)
noise2 = soundgen:::generateNoise(len = samplingRate * .15, noise = c(0, -80),
  rolloffNoise = c(4, -6), attackLen = 5)
noise3 = soundgen:::generateNoise(len = samplingRate * .25, noise = c(0, -40),
  rolloffNoise = c(4, -20), attackLen = 5)
# playme(c(noise2, noise3), samplingRate)

## Not run: 
playback = list(TRUE, FALSE, 'aplay', 'vlc')[[1]]
# 1.2 s of noise with rolloff changing from 0 to -12 dB above 2 kHz
noise = generateNoise(len = samplingRate * 1.2,
  rolloffNoise = c(0, -12), noiseFlatSpec = 2000,
  samplingRate = samplingRate, play = playback)
# spectrogram(noise, samplingRate)

# Similar, but using the dataframe format to specify a more complicated
# contour for rolloffNoise:
noise = generateNoise(len = samplingRate * 1.2,
  rolloffNoise = data.frame(time = c(0, .3, 1), value = c(-12, 0, -12)),
  noiseFlatSpec = 2000, samplingRate = samplingRate, play = playback)
# spectrogram(noise, samplingRate)

# To create a sibilant [s], specify a single strong, broad formant at ~7 kHz:
wl = 1024
formantFilter = getFormantFilter(
  nr = wl %/% 2 + 1, nc = 1, samplingRate = samplingRate,
 formants = list('f1' = data.frame(time = 0, freq = 7000,
                                   amp = 50, width = 2000)))
noise = fade(generateNoise(len = samplingRate,
  samplingRate = samplingRate, formantFilter = as.numeric(formantFilter),
  play = playback), samplingRate = samplingRate)
# plot(formantFilter, type = 'l')
meanSpectrum(noise, samplingRate)

# Low-frequency, wind-like noise
formantFilter = getFormantFilter(
  nr = 50, nc = 1, lipRad = 0,
  samplingRate = samplingRate, formants = list('f1' = list(
    freq = 250, amp = 30, width = 150)),
    formantDepStoch = 0, plot = TRUE)
noise = fade(generateNoise(len = samplingRate,
  samplingRate = samplingRate, formantFilter = as.numeric(formantFilter),
  play = playback))
spectrogram(noise, samplingRate, ylim = c(0, 2))

# Manual filter, e.g. for a kettle-like whistle (narrow-band noise)
formantFilter = c(rep(0, 100), 120, rep(0, 100))  # any length is fine
# plot(formantFilter, type = 'b')  # narrow-band filter at Nyquist / 2, here 4 kHz
noise = fade(generateNoise(len = samplingRate, formantFilter = formantFilter,
  samplingRate = samplingRate, play = playback))
spectrogram(noise, samplingRate)

# Compare to a similar sound created with soundgen()
# (aperiodic noise only, a single formant at 4 kHz)
noise_s = soundgen(pitch = NULL,
  noise = data.frame(time = c(0, 1000), value = c(0, 0)),
  formants = list(f1 = data.frame(freq = 4000, amp = 80, width = 20)),
  play = playback)

## End(Not run)

Get bandwidth

Description

Calculates formant bandwidth as a function of formant frequencies using a modified version of TMF-63 formula. Namely, above 500 Hz it follows the original formula from Tappert, Martony, and Fant (TMF)-1963, and below 500 Hz it applies a correction to allow for energy losses at low frequencies. See Khodai-Joopari & Clermont (2002), "Comparison of formulae for estimating formant bandwidths". Below 250 Hz the bandwidth is forced to drop again to avoid very large values near zero (just guesswork!)

Usage

getBandwidth(f)

Arguments

Examples

f = 1:5000
plot(f, soundgen:::getBandwidth(f), type = 'l',
  xlab = 'Formant frequency, Hz', ylab = 'Estimated bandwidth, Hz')

Get duration

Description

Returns the duration of one or more audio files (mostly useful for running on an entire folder). If threshold is set, it also removes the leading and trailing silences or near-silences, thus returning the duration of relatively loud central fragments of each sound. Silences are located based on the amplitude of root mean square (RMS) amplitude with getRMS. Note that the threshold is set relative to the observed maximum RMS, just as in analyze. This means that even very quiet sounds are not treated as pure silence.

Usage

getDuration(
  x,
  samplingRate = NULL,
  silence = 0.01,
  rms_pars = list(windowLength = 20, step = 5),
  reportEvery = NULL,
  cores = 1
)

Arguments

Value

For a single file, a one-row data frame with columns $duration and $duration_noSilence. For multiple files, a data frame with one row per file, including a $file column. If silence = NULL, duration_noSilence (duration without leading and trailing silence) is NA. If no RMS frame exceeds the threshold, duration_noSilence is 0.

See Also

analyze getLoudness

Examples

s = c(rep(0, 550), runif(400, -1, 1), rep(0, 50))
osc(s, samplingRate = 1000)
# true duration_noSilence is 400 ms
getDuration(s, samplingRate = 1000, silence = .01)
getDuration(s, samplingRate = 1000, silence = .1,
            rms_pars = list(windowLength = 5, step = 1))

## Not run: 
d = getDuration('~/Downloads/temp')
d
hist(d$duration - d$duration_noSilence)

## End(Not run)

Entropy

Description

Returns Wiener or Shannon entropy of an input vector such as the spectrum of a sound. Non-positive input values are converted to a small positive number (convertNonPositive). If all elements are zero, returns NA.

Usage

getEntropy(
  x,
  type = c("wiener", "shannon"),
  normalize = FALSE,
  convertNonPositive = 1e-10
)

Arguments

Value

The entropy of input vector as a single number.

Examples

# Here are four simplified power spectra, each with 9 frequency bins:
s = list(
  onePeak = c(rep(0, 4), 1, rep(0, 4)),         # a single peak in spectrum
  harmonics = c(0, 0, 1, 0, 0, .75, 0, 0, .5),  # perfectly periodic, with 3 harmonics
  silent = rep(0, 9),                           # a silent frame
  noise = rep(1, 9)                             # white noise
)

# Wiener entropy is ~0 for periodic, NA for silent, 1 for white noise
lapply(s, function(x) round(getEntropy(x), 2))

# Shannon entropy is ~0 for periodic with a single harmonic, moderate for
# periodic with multiple harmonics, NA for silent, highest for white noise
lapply(s, function(x) round(getEntropy(x, type = 'shannon'), 2))

# Normalized Shannon entropy - same but forced to be 0 to 1
lapply(s, function(x) round(getEntropy(x,
  type = 'shannon', normalize = TRUE), 2))

Get amplitude envelope

Description

Calculates a smoothed envelope of a waveform based on peaks, running average, root mean square (intensity), or envelope of the analytic signal with Hilbert transform. Algorithm: calculates one envelope value per frame of length wl, then upsamples to the original sampling rate with spline.

Usage

getEnv(
  x,
  method = c("rms", "hil", "peak", "mean"),
  wl = 200,
  overlap = 0,
  step = NULL,
  upsample = TRUE
)

Arguments

Value

The envelope as a numeric vector on the original scale. If upsample = TRUE, it has the same length as the input, regardless of the amount of smoothing.

Examples

a = rnorm(500) * seq(1, 0, length.out = 500)
wl = 50
scale = max(abs(a))
plot(a, type = 'l', ylim = c(-scale, scale))
lines(getEnv(a, 'rms', wl), col = 'red')
lines(getEnv(a, 'peak', wl), col = 'green')
lines(getEnv(a, 'hil', wl), col = 'blue')
lines(getEnv(a, 'mean', wl), lty = 3, lwd = 3)

# No upsampling (short output)
getEnv(1:16, 'mean', wl = 5, upsample = FALSE)
env_short = getEnv(a, 'rms', wl = wl, overlap = 50, upsample = FALSE)
plot(a, type = 'l')
lines(seq(1, length(a), length.out = length(env_short)), env_short, col = 'red')

Formant filter

Description

Prepares a frequency-domain filter for adding formants to a sound. Formants are specified as a list containing time, frequency, amplitude, and width values for each vocal tract resonance (see examples). For more information, see https://cogsci.se/soundgen/sound_generation.html.

Usage

getFormantFilter(
  nr,
  nc,
  formants = NA,
  formantDep = 1,
  formantWidth = 1,
  lipRad = 6,
  noseRad = 4,
  mouth = NA,
  mouthOpenThres = 0,
  openMouthBoost = 0,
  vocalTract = NULL,
  temperature = 0.025,
  formDrift = 0.3,
  formDisp = 0.2,
  formantDepStoch = 1,
  formantCeiling = NULL,
  samplingRate = 16000,
  speedSound = 35400,
  smoothing = list(interpol = "splineFC"),
  output = c("simple", "detailed"),
  plot = FALSE,
  duration = NULL,
  colorTheme = "bw",
  col = NULL,
  xlab = "Time",
  ylab = "Frequency, kHz",
  ...
)

Arguments

Value

Spectral filter on a linear scale (not dB): a matrix with frequency bins in rows and time steps in columns. Accordingly, rownames of the output give central frequency of each bin (in kHz), while colnames give time stamps (in ms if duration is specified, otherwise 0 to 1).

Examples

# [a] with only F1-F3 visible, with no stochasticity
e = getFormantFilter(nr = 512, nc = 50, duration = 300,
  formants = 'a', temperature = 0, plot = TRUE, col = heat.colors(150))
# image(t(e))  # to plot the output on a linear scale instead of dB

# some "wiggling" of specified formants plus extra formants on top
e = getFormantFilter(nr = 512, nc = 50,
  formants = c(860, 1430, 2900),
  temperature = 0.1, formantDepStoch = 1, plot = TRUE)

# a schwa based on variable length of vocal tract
e = getFormantFilter(nr = 512, nc = 50, formants = NA,
  vocalTract = list(time = c(0, .4, 1), value = c(13, 18, 17)),
  temperature = .1, plot = TRUE)

# no formants at all, only lip radiation
e = getFormantFilter(nr = 512, nc = 1, lipRad = 6,
  formants = NA, temperature = 0, plot = FALSE)
plot(e[, 1], type = 'l')              # linear scale
plot(20 * log10(e[, 1]), type = 'l')  # dB scale - 6 dB/oct

# mouth opening
e = getFormantFilter(nr = 512, nc = 50,
  vocalTract = 16, plot = TRUE, lipRad = 6, noseRad = 4,
  mouth = data.frame(time = c(0, .5, 1), value = c(0, 0, .5)))

# scale formant amplitude and/or bandwidth
e1 = getFormantFilter(nr = 512, nc = 1,
  formants = 'a', formantWidth = 1, formantDep = 1)  # defaults
e2 = getFormantFilter(nr = 512, nc = 1,
  formants = 'a', formantWidth = 1.5, formantDep = 1.5)
plot(as.numeric(rownames(e2)), 20 * log10(e2[, 1]),
     type = 'l', xlab = 'KHz', ylab = 'dB', col = 'red', lty = 2)
points(as.numeric(rownames(e1)), 20 * log10(e1[, 1]), type = 'l')

# manual specification of formants
e3 = getFormantFilter(
  nr = 512, nc = 50, samplingRate = 16000, plot = TRUE,
  formants = list(
    f1 = list(freq = c(900, 500), amp = c(30, 35), width = c(80, 50)),
    f2 = list(freq = c(1900, 2500), amp = c(25, 30), width = 100),
    f3 = list(freq = 3400, amp = 30, width = 120)
))

# extra zero-pole pair (doesn't affect estimated VTL and thus the extra
# formants added on top)
e4 = getFormantFilter(
  nr = 512, nc = 50, samplingRate = 16000, plot = TRUE,
  formants = list(
    f1 = list(freq = c(900, 500), amp = c(30, 35), width = c(80, 50)),
    f1.5 = list(freq = 1300, amp = -15),
    f1.7 = list(freq = 1500, amp = 15),
    f2 = list(freq = c(1900, 2500), amp = c(25, 30), width = 100),
    f3 = list(freq = 3400, amp = 30, width = 120)
))
plot(as.numeric(rownames(e4)), 20 * log10(e3[, ncol(e3)]),
     type = 'l', xlab = 'KHz', ylab = 'dB')
points(as.numeric(rownames(e4)), 20 * log10(e4[, ncol(e4)]),
       type = 'l', col = 'red', lty = 2)

Get HNR

Description

Calculates the harmonics-to-noise ratio (HNR), that is, the ratio of the intensity of the harmonic component to the intensity of the noise component, following Boersma (1993). Normally called internally by analyze, but can also be called directly on a time series.

Usage

getHNR(
  x = NULL,
  samplingRate = NA,
  acf_x = NULL,
  lag_min = 2,
  lag_max = NULL,
  interpol = c("sinc", "spline", "parab", "none"),
  wn = "hanning",
  idx_max = NULL,
  win_sinc = NULL
)

Arguments

Value

A list:

f0

frequency (Hz) corresponding to the ACF peak

max_acf

height of the ACF peak, 0 to 1

HNR

harmonics-to-noise ratio in dB: 10 * log10(max_acf / (1 - max_acf))

References

Boersma, P. (1993). Accurate short-term analysis of the fundamental frequency and the harmonics-to-noise ratio of a sampled sound. In Proceedings of the Institute of Phonetic Sciences (Vol. 17, No. 1193, pp. 97–110).

Examples

signal = sin(2 * pi * 150 * (1:16000) / 16000)
signal = signal / sqrt(mean(signal ^ 2))
noise = rnorm(16000)
noise = noise / sqrt(mean(noise ^ 2))
SNR = 40  # ground truth
s = signal + noise * 10 ^ (-SNR / 20)
soundgen:::getHNR(s, 16000, lag_min = 16000 / 1000,
  lag_max = 16000 / 75, interpol = 'none')
soundgen:::getHNR(s, 16000, lag_min = 16000 / 1000,
  lag_max = 16000 / 75, interpol = 'sinc')

Discrete random walk

Description

Takes a continuous random walk and converts it to continuous epochs of repeated values 0/1/2, each at least minLength points long. 0/1/2 correspond to different noise regimes: 0 = no noise, 1 = subharmonics, 2 = subharmonics and jitter/shimmer.

Usage

getIntegerRandomWalk(
  rw,
  nonlinBalance = 50,
  minLength = 50,
  q1 = NULL,
  q2 = NULL,
  plot = FALSE
)

Arguments

Value

A vector of integers (0/1/2) of the same length as rw.

Examples

rw = soundgen:::getRandomWalk(len = 100, rw_range = 100, rw_smoothing = .2)
r = soundgen:::getIntegerRandomWalk(rw, nonlinBalance = 75,
                         minLength = 10, plot = TRUE)
r = soundgen:::getIntegerRandomWalk(rw, nonlinBalance = 15,
                         q1 = 30, q2 = 70,
                         minLength = 10, plot = TRUE)

Get loudness

Description

Estimates subjective loudness and sharpness of audio. Based on EMBSD speech quality measure, particularly the MATLAB code in Yang (1999) and Timoney et al. (2004). Note that there are many ways to estimate loudness and many other factors, ignored by this model, that could influence subjectively experienced loudness. Please treat the output with a healthy dose of skepticism! Also note that the absolute value of calculated loudness critically depends on the chosen "measured" sound pressure level (SPL). getLoudness estimates how loud a sound will be experienced if it is played back at an SPL of SPL_measured dB. The most meaningful way to use the output is to compare the loudness of several sounds analyzed with identical settings or of different segments within the same recording.

Usage

getLoudness(
  x,
  samplingRate = NULL,
  scale = NULL,
  from = NULL,
  to = NULL,
  input = c("spec", "audSpec"),
  windowLength = 50,
  step = NULL,
  overlap = 50,
  SPL_measured = 70,
  spreadSpectrum = FALSE,
  sharpnessMethod = c("aures", "DIN45692", "bismarck"),
  summaryFun = c("mean", "median", "sd"),
  reportEvery = NULL,
  cores = 1,
  plot = TRUE,
  savePlots = FALSE,
  embed = FALSE,
  main = NULL,
  ylim = NULL,
  width = 900,
  height = 500,
  units = "px",
  res = NA,
  mar = c(5.1, 4.1, 4.1, 4.1),
  ...
)

Arguments

Details

Algorithm: calibrates the sound to the desired SPL (Timoney et al., 2004), extracts a spectrogram with frequencies on the bark scale, optionally spreads the spectrum to account for frequency masking across the critical bands (Yang, 1999), converts dB to phon by using standard equal loudness curves (ISO 226), converts phon to sone (Timoney et al., 2004), sums across all critical bands, and applies a correction coefficient to standardize output. Calibrated so as to return a loudness of 1 sone for a 1 kHz pure tone with SPL of 40 dB and spreadSpectrum = FALSE. Sharpness is calculated as the weighted first moment of specific loudness on the Bark scale.

Value

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

$detailed contains per-file results. If multiple sounds are analyzed, $detailed is a list of per-sound lists. If a single sound is analyzed, it is simplified to a single list. Each list contains:

loudness

a vector of loudness in sone units per STFT frame

specSone

spectrum in bark-sone: a matrix of loudness values in sone, with frequency on the bark scale in rows and time (STFT frames) in columns

loudnessPhon, specPhon

same in phon instead of sone units

sharpness

a vector of sharpness in acum units per STFT frame

audSpec

auditory spectrogram used to calculate loudness and sharpness

$summary is a dataframe of summary loudness measures (one row per file). If summaryFun is NULL, $summary is NULL.

References

• ISO 226 as implemented by Jeff Tackett (2005) on https://www.mathworks.com/matlabcentral/fileexchange/ 7028-iso-226-equal-loudness-level-contour-signal

• Timoney, J., Lysaght, T., Schoenwiesner, M., & MacManus, L. (2004). Implementing loudness models in matlab.

• Yang, W. (1999). Enhanced Modified Bark Spectral Distortion (EMBSD): An Objective Speech Quality Measure Based on Audible Distortion and Cognitive Model. Temple University.

See Also

getRMS analyze

Examples

sounds = list(
  noise_1KHz = soundgen:::zeroOne(bandpass(rnorm(8000), 16000,
    lwr = 900, upr = 1100)) * 2 -1, # narrow-band noise at 1 KHz
  white_noise = runif(8000, -1, 1),  # white noise
  white_noise2 = runif(8000, -1/2, 1/2),  # ~6 dB quieter
  pure_tone_1KHz = sin(2*pi*1000/16000*(1:8000)),  # pure tone at 1 kHz
  pure_tone_100Hz = sin(2*pi*100/16000*(1:8000))  # pure tone at 100 Hz
)
# playme(sounds)
l = getLoudness(
    x = sounds, samplingRate = 16000, scale = 1,
    windowLength = 40, step = NULL, input = c('spec', 'audSpec')[1],
    overlap = 50, SPL_measured = 60,
    plot = FALSE)
l$summary
# loudness depends on amplitude if "scale" is provided (cf. sounds 2 and 3)
# narrowband noise / tone at 1 kHz, 60 dB: sharpness ~=1 acum, loudness ~=4 sone

# a steady glissando from 125 to 8000 Hz (constant on a musical scale)
pitch = exp(seq(log(125), log(8000), length.out = 16000))
s = sinpi(2 * cumsum(pitch) / 16000)
l1 = getLoudness(s, samplingRate = 16000, SPL_measured = 70)
# steady SPL, but variable loudness

# The estimated loudness and sharpness depend on target SPL
l2 = getLoudness(s, samplingRate = 16000, SPL_measured = 40, plot = FALSE)
l1$summary$loudness_mean
l2$summary$loudness_mean

# ...but not (much) on windowLength and samplingRate
l3 = getLoudness(s, samplingRate = 16000, SPL_measured = 40,
  windowLength = 50, plot = FALSE)
l3$summary$loudness_mean

## Not run: 
# Using auditory spectrogram as input instead of STFT (slower)
l4 = getLoudness(s, samplingRate = 16000, SPL_measured = 40, input = 'audSpec')
l4$summary$loudness_mean

# Process all audio files in a folder
l5 = getLoudness('~/Downloads/temp', savePlots = TRUE)
l5$summary

## End(Not run)

Get peak frequency

Description

Performs STFT and finds the dominant frequency within a target range for each frame, together with its spectral purity (proportion of normalized spectral magnitude at the peak) and, when a meaningful DC component is present, the modulation depth derived from the AC/DC ratio. Used by "getAM_env" for measuring amplitude modulation and by analyze for measuring frequency modulation.

Usage

getPeakFreq(
  x,
  samplingRate,
  freqRange = NULL,
  overlap = 75,
  parab = TRUE,
  plot = FALSE
)

Arguments

Details

For a sinusoidal modulator with depth m = amDep / 100, the one-sided (un-doubled) envelope spectrum satisfies |X(f_am)| / |X(0)| = m / (4 - 2m), giving amDep = 400 * purity / (dc + 2 * purity). The sum-to-one normalization cancels in this ratio, so the formula is unaffected by spectral leakage.

Value

A dataframe with one row per STFT frame:

time

time stamp, ms

freq

peak frequency, Hz (NA for frames with no usable energy in the target frequency range)

purity

peak magnitude as a proportion of the frame's total spectral magnitude. Normally approximately 0 to 1, but not clamped; values outside this range are preserved as diagnostic information.

dep

approximate modulation depth, on a 0 to 100 scale for an ideal sinusoidal envelope. Not clamped; values outside 0-100 can occur for non-sinusoidal or otherwise uncalibrated inputs. NA if the DC component is zero or non-finite.

Examples

{
# White noise with sinusoidal AM
amFreq = 10; amDep = 60
am = .5 + .5 * cospi(pi * amFreq * (1:1000) / 4000)
env = rnorm(4000) * (1 - am * amDep / 100)
plot(env, type = 'l')

soundgen:::getPeakFreq(env, samplingRate = 4000, freqRange = c(5, 50))
}

Zero-crossing rate

Description

A less precise, but very quick method of pitch tracking based on measuring zero-crossing rate in bandpass-filtered audio. Recommended for processing long recordings with typical pitch values well below the first formant frequency, such as speech. Calling this function is considerably faster than using the same pitch-tracking method in analyze. Note that, unlike analyze(), it returns the times of individual zero crossings (hopefully corresponding to glottal cycles) instead of pitch values at fixed time intervals.

Usage

getPitchZc(
  x,
  samplingRate = NULL,
  scale = NULL,
  from = NULL,
  to = NULL,
  pitchFloor = 50,
  pitchCeiling = 400,
  zcThres = 0.1,
  zcWin = 5,
  silence = 0.04,
  envWin = 5,
  certMethod = c("variab", "autocor"),
  summaryFun = c("mean", "sd"),
  reportEvery = NULL
)

Arguments

Details

Algorithm: the audio is bandpass-filtered from pitchFloor to pitchCeiling, and the timing of all zero crossings is saved. This is not enough, however, because aperiodic sounds like white noise also have plenty of zero crossings. Accordingly, an attempt is made to detect voiced segments (or steady musical tones, etc.) by looking for stable regions, with several zero-crossings at relatively regular intervals (see parameters zcThres and zcWin). Very quiet parts of audio are also treated as not having a pitch.

Value

A list with descriptives per file (@summary) and per frame (@detailed), including

time

time stamps of all zero crossings except the last one, after bandpass-filtering

pitch

pitch calculated from the time between consecutive zero crossings

cert

certainty in each pitch candidate calculated from local pitch stability, 0 to 1

See Also

analyze

Examples

data(speechEx, package = 'soundgen')
# spectrogram(speechEx)
zc = soundgen:::getPitchZc(speechEx, pitchCeiling = 250)
plot(zc$detailed[, c('time', 'pitch')], type = 'b')
spectrogram(speechEx, extraContour = zc$detailed$pitch, ylim = c(0, 2))

## Not run: 
# process all files in a folder
zc = soundgen:::getPitchZc('~/Downloads/temp')
zc$summary

## End(Not run)

RMS amplitude

Description

Calculates root mean square (RMS) amplitude in overlapping windows, providing an envelope of sound intensity. Longer windows provide smoother, more robust estimates; shorter windows and more overlap improve temporal resolution, but they also increase processing time and make the contour less smooth.

Usage

getRMS(
  x,
  samplingRate = NULL,
  scale = NULL,
  from = NULL,
  to = NULL,
  windowLength = 50,
  step = NULL,
  overlap = 70,
  stereo = c("left", "right", "average", "both"),
  killDC = FALSE,
  normalize = TRUE,
  windowDC = 200,
  summaryFun = "mean",
  reportEvery = NULL,
  cores = 1,
  plot = FALSE,
  savePlots = FALSE,
  embed = FALSE,
  main = NULL,
  xlab = "",
  ylab = "",
  type = "b",
  col = "green",
  lwd = 2,
  width = 900,
  height = 500,
  units = "px",
  res = NA,
  ...
)

Arguments

Details

Note that you can also get similar estimates per frame from analyze on a normalized scale of 0 to 1, but getRMS is much faster, operates on the original scale, and plots the amplitude contour. If you need RMS for the entire sound instead of per frame, you can simply calculate it as sqrt(mean(x^2)), where x is your waveform. Having RMS estimates per frame gives more flexibility: RMS per sound can be calculated as the mean / median / max of RMS values per frame.

Value

A list containing:

$detailed:

a list of RMS amplitudes per frame for each sound, on the scale of input; names give time stamps for the center of each frame, in ms.

$summary:

a dataframe with summary measures, one row per sound

See Also

analyze getLoudness

Examples

s = soundgen() + .25  # with added DC offset
# osc(s)
r = getRMS(s, samplingRate = 16000, from = .05,
  windowLength = 40, overlap = 50, killDC = TRUE,
  plot = TRUE, type = 'l', lty = 2, main = 'RMS envelope')
r

# short window = jagged envelope
r = getRMS(s, samplingRate = 16000,
  windowLength = 5, overlap = 0, killDC = TRUE,
  plot = TRUE, col = 'blue', pch = 13, main = 'RMS envelope')

 # stereo
 wave_stereo = tuneR::Wave(
   left = runif(1000, -1, 1) * 16000,
   right = runif(1000, -1, 1) / 3 * 16000,
   bit = 16, samp.rate = 4000)
 getRMS(wave_stereo)$summary
 getRMS(wave_stereo, stereo = 'right')$summary
 getRMS(wave_stereo, stereo = 'average')$summary
 getRMS(wave_stereo, from = .05,
   stereo = 'both', plot = TRUE)$summary

## Not run: 
r = getRMS('~/Downloads/temp', savePlots = TRUE)
r$summary

# Compare:
analyze('~/Downloads/temp', pitchMethods = NULL,
        plot = FALSE)$summary$ampl_mean
# (per STFT frame, but should be very similar)

# User-defined summary functions:
ran = function(x) diff(range(x))
meanSD = function(x) {
  paste0('mean = ', round(mean(x), 2), '; sd = ', round(sd(x), 2))
}
getRMS('~/Downloads/temp', summaryFun = c('mean', 'ran', 'meanSD'))$summary

## End(Not run)

Random walk

Description

Generates a random walk with flexible control over its range, trend, and smoothness. It works by calling stats::rnorm at each step and taking a cumulative sum of the generated values. Smoothness is controlled by initially generating a shorter random walk and upsampling.

Usage

getRandomWalk(
  len,
  rw_range = 1,
  rw_smoothing = 0.2,
  method = c("linear", "spline")[2],
  trend = 0
)

Arguments

Value

A numeric vector of length len and range from 0 to rw_range.

Examples

plot(soundgen:::getRandomWalk(len = 1000, rw_range = 5, rw_smoothing = 0))
plot(soundgen:::getRandomWalk(len = 1000, rw_range = 5, rw_smoothing = .2))
plot(soundgen:::getRandomWalk(len = 1000, rw_range = 5, rw_smoothing = .95))
plot(soundgen:::getRandomWalk(len = 1000, rw_range = 5, rw_smoothing = .99))
plot(soundgen:::getRandomWalk(len = 1000, rw_range = 5, rw_smoothing = 1))
plot(soundgen:::getRandomWalk(len = 1000, rw_range = 15,
  rw_smoothing = .2, trend = c(.1, -.1)))
plot(soundgen:::getRandomWalk(len = 1000, rw_range = 15,
  rw_smoothing = .2, trend = c(15, -1)))

Rolloff of harmonics in source spectrum

Description

A periodic sound often contains many harmonics of the fundamental frequency (f0). For instance, soundgen generates voice source as a combination of pure tones - one for each partial in the source spectrum. Harmonics typically become weaker at higher frequencies, and the rolloff (roll-off) function describes the loss of energy in upper harmonics relative to f0. getRolloff provides flexible control over this rolloff function: simple exponential (rolloff), rolloffOct that varies per octave, and rolloffKHz for rolloff correction depending on f0. It calculates the rolloff function based on these parameters and returns the amplitude of each harmonic.

Usage

getRolloff(
  pitch_per_gc = c(440),
  nHarmonics = 1000,
  rolloff = -6,
  rolloffOct = 0,
  rolloffKHz = 0,
  baseline = 200,
  dynamicRange = 80,
  samplingRate = 16000,
  plot = FALSE
)

Arguments

Value

A matrix of amplitude multiplication factors for adjusting the amplitude of harmonics relative to f0 (1 = no adjustment, 0 = silent). Each row of output contains one harmonic, and each column contains one glottal cycle.

See Also

soundgen

Examples

# steady exponential rolloff of -12 dB per octave
rolloff1 = getRolloff(pitch_per_gc = 150, rolloff = -12,
  rolloffOct = 0, rolloffKHz = 0, plot = TRUE)
# the rate of rolloff slows down by 1 dB each octave
r = getRolloff(pitch_per_gc = 100, rolloff = -6, rolloffOct = 0)
r1 = getRolloff(pitch_per_gc = 100, rolloff = -6, rolloffOct = -1)
round(20*log10(r[c(1, 2, 4, 8, 16)]), 1)
round(20*log10(r1[c(1, 2, 4, 8, 16)]), 1)

# rolloff can be made to depend on f0 using rolloffKHz
rolloff = getRolloff(pitch_per_gc = c(150, 400, 800),
  rolloffOct = 0, rolloffKHz = -3, plot = TRUE)
# without the correction for f0 (rolloffKHz),
  # high-pitched sounds have the same rolloff as low-pitched sounds,
  # producing unnaturally strong high-frequency harmonics
rolloff = getRolloff(pitch_per_gc = c(150, 400, 800),
  rolloffOct = 0, rolloffKHz = 0, plot = TRUE)

# dynamic rolloff (varies over time)
rolloff = getRolloff(pitch_per_gc = c(150, 250),
                     rolloff = c(-12, -18, -24), plot = TRUE)

Calculate prop or fluctuation

Description

Calculates psychoacoustic roughness or slow fluctuation strength as the proportion of a modulation spectrum that falls within a specified range of temporal modulation frequencies - both as a single total value and per frequency band. Roughness or fluctuation can be calculated either as the simple proportion of the spectrum within modRange, or as a weighted version using a lognormal weighting function centered at modMean with spread modSD. Recommended weighting function for measuring roughness: modMean = 100, modSD = 8 (100 Hz ± 8 semitones); for measuring low-frequency fluctuation: modMean = 4, modSD = 12 (4 Hz ± 12 semitones). Because the result is a proportion (numerator and denominator are summed over the same modulation-frequency columns), positive and negative temporal modulation frequencies do not need to be folded: any mirroring cancels out.

Usage

getRoughness(m, modRange = NULL, modMean = 100, modSD = 8, plot = FALSE)

Arguments

Value

A list with two components:

• prop: the proportion of MS in the target range of modulation frequencies, in percent (NA if no valid modRange or weighting parameters, empty input, or no energy in the spectrum).

• details: a data.frame with one row per frequency band and columns freq (band label), prop (band-wise contribution to total proportion, %), and amp (mean band amplitude, normalized to the maximum band).

See Also

modulationSpectrum

Examples

s = soundgen(jitterDep = 2, amFreq = 80, amDep = 70, addSilence = 0)
m = modulationSpectrum(s, 16000, amRes = NULL, plot = FALSE)$detailed$original

## Ex. 1: simple proportion within modRange
r = soundgen:::getRoughness(m, modRange = c(30, Inf), plot = TRUE)
r = soundgen:::getRoughness(m, modRange = c(30, 150), plot = TRUE)
plot(r$details$freq, r$details$prop, type = 'b')
r$prop == sum(r$details$prop)

## Ex. 2: lognormal weighting function instead of modRange
# narrow
soundgen:::getRoughness(m, modRange = NULL, modMean = 75,
  modSD = 3, plot = TRUE)
# broad
soundgen:::getRoughness(m, modRange = NULL, modMean = 75,
  modSD = 10, plot = TRUE)
# truncated at modRange
soundgen:::getRoughness(m, modRange = c(30, 200), modMean = 75,
  modSD = 10, plot = TRUE)
# NB: a very large modSD means we essentially revert to a simple proportion

## Ex. 3: measure low-frequency AM fluctuation instead of prop
# (just use a different weighting function)
s2 = sinpi(2*440*(1:16000)/16000) * (1 - cospi(2*5*(1:16000)/16000))
osc(s2, 16000)
m2 = modulationSpectrum(s2, 16000, amRes = NULL, plot = FALSE)$detailed$original
fl = soundgen:::getRoughness(m2, modRange = c(0.25, 30), modMean = 4,
  modSD = 12, plot = TRUE)
fl$prop  # in this case, a measure of slow fluctuation

Get sharpness

Description

Takes a sone-spectrogram and calculates the sharpness of each frame as the first moment of specific loudness on the Bark scale:

S = 0.11 * \frac{\int z * g(z) * N'(z) dz}{\int N'(z) dz}

By default, it is assumed that each row corresponds to 1 Bark. If rows are spaced differently (e.g. an auditory filter bank), supply bark and barkWidth.

Usage

getSharpness(
  specSone,
  method = c("aures", "DIN45692", "bismarck"),
  loudness = NULL,
  bark = NULL,
  barkWidth = NULL
)

Arguments

Value

A numeric vector of sharpness values in acum, one per STFT frame.

Examples

s1 = cos(2*pi*1000*(1:8000)/16000)
l1 = getLoudness(s1, 16000, SPL_measured = 60, sharpnessMethod = 'aures')
# image(t(l1$detailed$specSone))
l1$detailed$sharpness
soundgen:::getSharpness(l1$detailed$specSone, method = 'aures') # same

s2 = soundgen()
l2 = getLoudness(s2, 16000, SPL_measured = 60, sharpnessMethod = 'aures')
# image(t(l2$detailed$specSone))
plot(l2$detailed$sharpness, type = 'b', ylim = c(0, 5))
points(soundgen:::getSharpness(l2$detailed$specSone, method = 'DIN45692'),
  type = 'l', col = 'blue')
points(soundgen:::getSharpness(l2$detailed$specSone, method = 'bismarck'),
  type = 'l', col = 'green')
points(soundgen:::getSharpness(l2$detailed$specSone, method = 'aures'),
  type = 'l', col = 'orange')

Smooth contour from anchors

Description

Returns a smooth contour based on an arbitrary number of anchors - points through which we wish to draw a curve, and which may or may not be equidistant along the x (time) dimension. getSmoothContour is used extensively by soundgen for generating intonation contour, mouth opening, etc. Pitch contours are treated as a special case: values are log-transformed prior to smoothing, so that with 2 anchors we get a linear transition on a log scale (as if we were operating with musical notes rather than frequencies in Hz).

Usage

getSmoothContour(
  anchors = data.frame(time = c(0, 1), value = c(0, 1)),
  len = NULL,
  thisIsPitch = FALSE,
  normalizeTime = TRUE,
  interpol = "splineFC",
  discontThres = 0.05,
  jumpThres = 0.01,
  valueFloor = NULL,
  valueCeiling = NULL,
  plot = FALSE,
  xlim = NULL,
  ylim = NULL,
  xlab = "Time, ms",
  ylab = ifelse(thisIsPitch, "Frequency, Hz", "Amplitude"),
  main = ifelse(thisIsPitch, "Pitch contour", ""),
  samplingRate = 16000,
  voiced = NULL,
  contourLabel = NULL,
  NA_to_zero = TRUE,
  ...
)

Arguments

Value

A numeric vector of length len.

Examples

# long format: anchors are a dataframe
a = soundgen:::getSmoothContour(anchors = data.frame(
  time = c(50, 137, 300), value = c(0.03, 0.78, 0.5)),
  normalizeTime = FALSE,
  voiced = 200, valueFloor = 0, plot = TRUE, main = '',
  samplingRate = 16000) # breathing

# short format: anchors are a vector (equal time steps assumed)
a = soundgen:::getSmoothContour(anchors = c(350, 800, 600),
  len = 5500, thisIsPitch = TRUE, plot = TRUE,
  samplingRate = 3500) # pitch

# a single anchor gives constant value
a = soundgen:::getSmoothContour(anchors = 800,
  len = 500, thisIsPitch = TRUE, plot = TRUE, samplingRate = 500)

# two pitch anchors give loglinear F0 change
a = soundgen:::getSmoothContour(anchors = c(220, 440),
  len = 500, thisIsPitch = TRUE, plot = TRUE, samplingRate = 500)

## Two closely spaced anchors produce a pitch jump
# one curve through all anchors
a1 = soundgen:::getSmoothContour(anchors = list(time = c(0, .15, .2, .7, 1),
    value = c(360, 116, 550, 700, 610)), len = 500, thisIsPitch = TRUE,
    plot = TRUE, samplingRate = 500)
# two segments with a linear transition
a2 = soundgen:::getSmoothContour(anchors = list(time = c(0, .15, .17, .7, 1),
    value = c(360, 116, 550, 700, 610)), len = 500, thisIsPitch = TRUE,
    plot = TRUE, samplingRate = 500)
# two segments with an abrupt jump
a3 = soundgen:::getSmoothContour(anchors = list(time = c(0, .15, .155, .7, 1),
    value = c(360, 116, 550, 700, 610)), len = 500, thisIsPitch = TRUE,
    plot = TRUE, samplingRate = 500)
# compare:
plot(a2)
plot(a3)  # NB: the segment before the jump is upsampled to compensate

## Control the amount of smoothing
# linear - no smoothing
soundgen:::getSmoothContour(c(1, 3, 9, 10, 9, 9, 2), len = 100, plot = TRUE,
  interpol = 'linear')
# spline (smooth, but may overshoot slightly)
soundgen:::getSmoothContour(c(1, 3, 9, 10, 9, 9, 2), len = 100, plot = TRUE,
  interpol = 'splineFC')
# cosine - goes strictly through anchors, no overshoot, but possible ripples
soundgen:::getSmoothContour(c(1, 3, 9, 10, 9, 9, 2), len = 100, plot = TRUE,
  interpol = 'cosine')
# see ?interpolate for more options

## Upsample preserving leading and trailing NAs
anchors = data.frame(time =  c(1,  4,  5,  7,  10, 20, 23, 25, 30),
                     value = c(NA, NA, 10, 15, 12, NA, 17, 15, NA))
plot(anchors, type = 'b')
anchors_ups = soundgen:::getSmoothContour(
  anchors, len = 200,
  interpol = 'linear',
  NA_to_zero = FALSE,   # preserve NAs
  discontThres = 0)     # don't break into sub-contours
plot(anchors_ups, type = 'b')

Universal spectrogram extractor

Description

Routes user-friendly strings to various functions for producing spectrogram-like representations (spectrogram, mel-spectrograms, RMS envelopes, spectra, etc.), executes them, and standardizes the output to a matrix with features in rows and time (in ms) in columns.

Usage

getSpec(audio, specFun = "melfcc", specFun_pars = list())

Arguments

Value

Returns a spectrogram-like matrix with time in columns (ms) and features in rows (names as passed by specFun).

Examples

a = rnorm(100)
# call function by name as a string (must be among the recognized strings)
soundgen:::getSpec(audio = soundgen:::readAudio(a, samplingRate = 1600),
        specFun = 'powspec',  # adds freq labels
        specFun_pars = list(wintime = .02, steptime = .01))
soundgen:::getSpec(audio = soundgen:::readAudio(a, samplingRate = 1600),
        specFun = 'spectrogram',
        specFun_pars = list(windowLength = 25))
soundgen:::getSpec(audio = soundgen:::readAudio(a, samplingRate = 1600),
        specFun = 'stft',
        specFun_pars = list(wl = 20))
soundgen:::getSpec(audio = soundgen:::readAudio(a, samplingRate = 1600),
        specFun = 'melspec',
        specFun_pars = list(nbands = 5))
soundgen:::getSpec(audio = soundgen:::readAudio(a, samplingRate = 1600),
        specFun = 'env',
        specFun_pars = list(wl = 20))
soundgen:::getSpec(audio = soundgen:::readAudio(a, samplingRate = 1600),
        specFun = 'ssm',  # ssm itself calls getSpec()
        specFun_pars = list(specFun = 'stft', specFun_pars = list(wl = 20)))
soundgen:::getSpec(audio = soundgen:::readAudio(a, samplingRate = 1600),
        specFun = 'env',
        specFun_pars = list(wl = 20))
soundgen:::getSpec(audio = soundgen:::readAudio(a, samplingRate = 1600),
        specFun = 'env',
        specFun_pars = list(wl = 20))
soundgen:::getSpec(audio = soundgen:::readAudio(a, samplingRate = 1600),
        specFun = 'env',
        specFun_pars = list(wl = 20))
soundgen:::getSpec(audio = soundgen:::readAudio(a, samplingRate = 1600),
        specFun = 'spectrum',
        specFun_pars = list(wl = 20))

# pass a function object directly (must accept audio as a numeric vector)
soundgen:::getSpec(audio = soundgen:::readAudio(a, samplingRate = 1600),
        specFun = tuneR::powspec,
        specFun_pars = list(sr = 1600, wintime = .02, steptime = .01))
# (NB: frequency labels not added - must be handled by the specFun itself)
soundgen:::getSpec(audio = soundgen:::readAudio(a, samplingRate = 1600),
        specFun = spectrogram,
        specFun_pars = list(samplingRate = 1600, windowLength = 2, plot = FALSE))

# construct a custom function
soundgen:::getSpec(audio = soundgen:::readAudio(1:10, samplingRate = 20),
        specFun = function(x) matrix(cumsum(x), nrow = 1),
        specFun_pars = list())

Get spectral envelope

Description

Calculates a smoothed envelope of a magnitude spectrum or of each column of a spectrogram. This is good for removing the fine structure produced by harmonics of f0 and leaving only the overall spectral contour produced by resonances (formants). This is the source-filter separation step used by shiftFormants. All methods except "peak" smooth the log-magnitude spectrum and return the envelope on the original linear scale, with the same dimensions as the input.

Usage

getSpecEnv(
  spec,
  freqs = NULL,
  freqWindow = NULL,
  freqWindow_bins = NULL,
  method = c("cepstral", "gauss", "movavg", "peak"),
  plot = FALSE
)

Arguments

Details

The amount of smoothing is controlled by freqWindow_bins, which should normally equal the expected spacing between harmonics (f0, in bins): spectral details that vary on a faster scale are treated as source fine structure and removed, while slower variations (formants) are retained. For high-pitched or variable calls, increase freqWindow_bins to smooth more. Methods:

• "cepstral" (default): Gaussian low-pass liftering of the real cepstrum (FFT of the log spectrum). The harmonic ripple, which has a period of freqWindow_bins bins, is attenuated to ~exp(-4) = 2 amplitude, while the formant envelope is preserved. By the convolution theorem, this is equivalent to Gaussian smoothing of the log spectrum with SD = freqWindow_bins / 2, but computed in the quefrency domain.

• "gauss": Gaussian blur of the log spectrum along the frequency axis with SD = freqWindow_bins / 2 - the same low-pass filter as "cepstral", but with explicit edge padding instead of circular (wrap-around) filtering.

• "movavg": moving average of the log spectrum with a rectangular window of width freqWindow_bins (rounded up to an odd number).

• "peak": moving maximum (morphological dilation), i.e. the upper envelope of the spectrum with a window of width freqWindow_bins. Much slower than the other methods on long inputs.

The edges of the spectrum are handled by padding (repeating the edge values) before smoothing and trimming afterwards, so the output always has the same length as the input. If spec is a matrix (rows = frequency bins, columns = time frames), all frames are smoothed at once.

Value

The spectral envelope on the original (linear magnitude) scale, as a numeric vector or matrix with the same dimensions as spec.

See Also

getEnv for the temporal envelope of a waveform; shiftFormants, which uses getSpecEnv to shift formants

Examples

# Synthetic spectrum: three formants plus harmonics 20 bins apart
N = 512
freq = 1:N
true_envelope = exp(-.5 * ((freq - 100) / 20)^2) +
  exp(-.5 * ((freq - 250) / 30)^2) +
  exp(-.5 * ((freq - 400) / 40)^2)
spectrum = true_envelope * (0.5 + 0.5 * abs(sin(pi * freq / 20)))

plot(freq, spectrum, type = 'l', log = 'y',
     main = 'Spectral envelope', xlab = 'Frequency, bins')
lines(freq, true_envelope, col = 'grey60', lwd = 4)
lines(getSpecEnv(spectrum, freqWindow_bins = 20, method = 'cepstral'),
      col = 'red', lwd = 2)
lines(getSpecEnv(spectrum, freqWindow_bins = 20, method = 'gauss'),
      col = 'orange', lwd = 2, lty = 2)
lines(getSpecEnv(spectrum, freqWindow_bins = 20, method = 'movavg'),
      col = 'blue', lwd = 2, lty = 3)
lines(getSpecEnv(spectrum, freqWindow_bins = 20, method = 'peak'),
      col = 'green', lwd = 2, lty = 4)
legend('bottom',
       legend = c('raw', 'truth', 'cepstral', 'gauss', 'movavg', 'peak'),
       col = c('black', 'grey60', 'red', 'orange', 'blue', 'green'),
       lwd = c(1, 4, 2, 2, 2, 2), lty = c(1, 1, 1, 2, 3, 4), bty = 'n')

# Smoothed spectral envelope of a single vowel
data(speechEx, package = 'soundgen')
spec = spectrum(speechEx, from = .15, to = .3, plot = FALSE)
env = getSpecEnv(spec, freqWindow = 500, plot = TRUE)

# Smooth a whole spectrogram at once (matrix input):
spec = stft_simple(speechEx@left[1:16000],
                   samplingRate = speechEx@samp.rate,
                   wl = 512, step = 256)
spec = Mod(spec[1:(nrow(spec) %/% 2 + 1), ])
env = getSpecEnv(spec, freqWindow = 500, plot = TRUE)

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

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)

Hilbert transforms

Description

hilbert_exact treats the input as one period of an infinite periodic signal (circular convolution), while hilbert_approx pads the input on both sides to a good length for FFT, which is faster at the cost of possible slight artifacts at the edges (~linear convolution). If the input is <3 samples long, the envelope is calculated simply as Mod(x). Note: only hilbert_approx is fast enough to be used for extracting the envelopes of typical audio.

Usage

hilbert_exact(x)

hilbert_approx(x)

Arguments

Value

A list with two components:

hilbert

Hilbert transform

envelope

Instantaneous amplitude of the analytic signal (envelope)

Examples

# signal: amplitude-modulated sine wave
t = seq(0, 1, length.out = 477)
carrier = cos(2 * pi * 50 * t)
modulator = 1 + 0.5 * cos(2 * pi * 5 * t)
s = modulator * carrier  # amplitude modulated signal
hil_exact = hilbert_exact(s)
hil_approx = hilbert_approx(s)
plot(s, type = 'l')
points(hil_exact$envelope, type = 'l', col = 'blue')
points(hil_approx$envelope, type = 'l', col = 'red')

Formants in American vowels

Description

Typical relative frequencies of the first four formants measured in dF units (average spacing between formants, or formant dispersion) above or below schwa based on estimated VTL in American English, from Hillenbrand (1995), who measured F1-F4 in ~1.5K recordings (139 speakers, 12 vowels from each). Audio and formant measurements are freely available online: https://homepages.wmich.edu/~hillenbr/voweldata.html. The dataset below is the result of modeling Hillenbrand's data with brms: mvbind(F1rel, F2rel) ~ vowel + (vowel|speaker). It shows the most credible location of each vowel centroid in the F1Rel-F2Rel space.

Usage

hillenbrand

Format

An object of class data.frame with 12 rows and 5 columns.

Details

A dataframe of 12 observations and 5 columns: "vowel" = vowel (American English), "F1Rel" to "F4Rel" = formant frequencies in dF relative to their neutral, equidistant positions in a perfectly cylindrical vocal tract. See schwa - this is what schwa() returns as $ff_relative_dF

References

Hillenbrand, J., Getty, L. A., Clark, M. J., & Wheeler, K. (1995). Acoustic characteristics of American English vowels. The Journal of the Acoustical society of America, 97(5), 3099-3111.

Examples

plot(hillenbrand$F1Rel, hillenbrand$F2Rel, type = 'n')
text(hillenbrand$F1Rel, hillenbrand$F2Rel, labels = hillenbrand$vowel)

Interpolate matrix

Description

Performs a chosen type of separable interpolation across both rows and columns of a matrix, in effect up- or downsampling a matrix to required dimensions. Rownames and colnames are also interpolated as needed. Make sure there are no NAs in the input to avoid unpredictable behavior (dropped by approx and spline, propagated by resample).

Usage

interpolMatrix(
  m,
  nr = NULL,
  nc = NULL,
  interpol = c("approx", "spline", "resample")
)

Arguments

Examples

m = matrix(1:12 + rnorm(12, 0, .2), nrow = 3)
rownames(m) = 1:3; colnames(m) = 1:4
soundgen:::interpolMatrix(m)  # just returns the original
soundgen:::interpolMatrix(m, nr = 10, nc = 7)
soundgen:::interpolMatrix(m, nr = 10, nc = 7, interpol = 'spline')
soundgen:::interpolMatrix(m, nr = 2, nc = 7)
soundgen:::interpolMatrix(m, nr = 2, nc = 3)

# input matrices can have a single row/column
soundgen:::interpolMatrix(matrix(1:5, nrow = 1), nc = 9)
soundgen:::interpolMatrix(matrix(1:5, ncol = 1), nr = 5, nc = 3)

# use resample() with a low-pass filter to avoid artifacts
a = matrix(c(rep(0, 10), 1, rep(0, 10)))
soundgen:::interpolMatrix(a, nr = 5)  # ok
soundgen:::interpolMatrix(a, nr = 6)  # spike lost
soundgen:::interpolMatrix(a, nr = 6, interpol = 'resample')  # ok

# complex values
cm = matrix(complex(real = 1:9, imaginary = 1:9), nrow = 3)
soundgen:::interpolMatrix(cm, nr = 5)

Interpolation

Description

Interpolates unevenly spaced points into a relatively smooth curve. If the points are evenly spaced but there are missing values and/or aliasing should be avoided, use resample instead.

Usage

interpolate(
  x,
  y,
  xout,
  method = c("splineFC", "spline", "constant", "linear", "approxLowPass", "sgolay",
    "pchip", "cosine", "cardinal", "hermite", "loess"),
  plot = FALSE,
  ...
)

interpol_approxLowPass(bandwidth = 0.1)

interpol_sgolay(p = 3, n = 15)

interpol_cardinal(tension = 0.5)

interpol_loess(span = 0.5)

Arguments

Value

A numeric vector of interpolated y‑values at the requested xout locations.

Interpolation methods

constant

Constant interpolation via approx. Fast, but no smoothing.

linear aka approx

Linear interpolation via approx. Fast, but not smooth.

spline

Cubic spline interpolation (FMM method) via spline. Fast, but overshoots.

splineFC

Monotone cubic interpolation using the Fritsch‑Carlson method (see splinefun). Moderately fast, less overshooting than the FMM spline.

approxLowPass

Linear interpolation followed by low‑pass filtering. Fast, smooth, but reduced range compared to original y. Constructor: interpol_approxLowPass(bandwidth).

sgolay

Linear interpolation followed by Savitzky‑Golay smoothing (see sgolayfilt). Fairly similar to approxLowPass, but much slower. Constructor: interpol_sgolay(p, n).

pchip

Piecewise Cubic Hermite Interpolating Polynomial (preserves monotonicity). Calls interp1 with method = "pchip".

cosine

Cosine‑eased (smoothstep) interpolation. Eases between anchor points using a cosine curve. Fast, passes exactly through input points.

cardinal

Cardinal spline interpolation (a generalization of Catmull‑Rom). Fast, can be forced to pass exactly through input points. Constructor: interpol_cardinal(tension).

hermite

Hermite spline interpolation with forced zero slope at local extrema to prevent overshoot. Fast, passes exactly through input points.

loess

Locally estimated scatterplot smoothing (LOESS) via loess. Smooth, but slow; may overshoot. Constructor: interpol_loess(span).

See Also

interpolateNA resample

Examples

x  = c(0, .15, .2, .3, .7, 1)
y  = c(360, 116, 550, 350, 700, 610)
xout = seq(0, 1, length.out = 100)

# Compare inbuilt interpolation methods
ms = c('constant', 'linear', 'spline', 'splineFC', 'approxLowPass',
        'sgolay', 'cosine', 'cardinal', 'hermite', 'loess')
op = par(c('mfrow', 'mar')); par(mfrow = c(4, 3), mar = c(2, 2, 3, 1))
for (m in ms) {interpolate(x, y, xout, method = m, plot = TRUE); title(m)}
par(op)

# Passing method‑specific parameters via ...
interpolate(x, y, xout, method = 'cardinal', tension = 0, plot = TRUE)
interpolate(x, y, xout, method = 'loess', span = 0.2, plot = TRUE)
interpolate(x, y, xout, method = 'loess', span = 0.9, plot = TRUE)

# Equivalent: passing a constructor (useful when forwarding through
# higher‑level functions like resample())
interpolate(x, y, xout, method = interpol_cardinal(tension = 0), plot = TRUE)
interpolate(x, y, xout, method = interpol_loess(span = 0.2), plot = TRUE)

# Passing a fully custom function
interpolate(x, y, xout, plot = TRUE,
  method = function(...) spline(x, y, xout = xout, method = 'natural')$y)

Interpolate NAs

Description

Takes a numeric vector, such as a pitch contour, and fills in the NAs, first by linear interpolation in the middle and then by constant or linear interpolation at the ends.

Usage

interpolateNA(x, idx_na = NULL, nPoints = 1)

Arguments

Value

The input numeric vector with NAs filled in by interpolation.

Examples

a = c(NA, 405, NA, 460, NA, NA, NA, 480, 490, NA, NA)
interpolateNA(a)
interpolateNA(a, nPoints = 3)
plot(interpolateNA(a), type = 'l', col = 'blue'); points(a)
plot(interpolateNA(a, nPoints = 3), type = 'l', col = 'blue'); points(a)
# Compare to approx - terminal NAs are simply trimmed
approx(a, na.rm = TRUE, n = length(a))$y

Invert spectrogram

Description

Transforms a spectrogram into a time series with inverse STFT. The problem is that an ordinary spectrogram preserves only the magnitude (modulus) of the complex STFT, while the phase is lost, and without phase it is impossible to reconstruct the original audio accurately. So there are a number of algorithms for "guessing" the phase that would produce an audio whose magnitude spectrogram is very similar to the target spectrogram. Useful for certain filtering operations that modify the magnitude spectrogram followed by inverse STFT, such as filtering in the spectro-temporal modulation domain.

Usage

invertSpectrogram(
  spec,
  samplingRate = NULL,
  windowLength,
  step,
  overlap = NULL,
  wn,
  specScale = c("spec", "power", "log", "dB"),
  initialPhase = c("spsi", "random", "zero"),
  nIter = 50,
  normalize = TRUE,
  play = FALSE,
  verbose = FALSE,
  plotError = TRUE
)

Arguments

Details

Algorithm: takes the spectrogram, makes an initial guess at the phase (zero, noise, or a more intelligent estimate by the SPSI algorithm), fine-tunes over nIter iterations with the GL algorithm, reconstructs the complex spectrogram using the best phase estimate, and performs inverse STFT. The single-pass spectrogram inversion (SPSI) algorithm is implemented as described in Beauregard et al. (2015) following the python code at https://github.com/lonce/SPSI_Python. The Griffin-Lim (GL) algorithm is based on Griffin & Lim (1984).

Value

Reconstructed audio as a numeric vector.

References

• Griffin, D., & Lim, J. (1984). Signal estimation from modified short-time Fourier transform. IEEE Transactions on Acoustics, Speech, and Signal Processing, 32(2), 236-243.

• Beauregard, G. T., Harish, M., & Wyse, L. (2015, July). Single pass spectrogram inversion. In 2015 IEEE International Conference on Digital Signal Processing (DSP) (pp. 427-431). IEEE.

See Also

spectrogram filterSoundByMS

Examples

# Create a spectrogram. NB: do NOT use zero-padding
samplingRate = 16000
windowLength = 40
step = 5
wn = 'gaussian'
# NB: the more detailed a spectrogram, the more precisely it can be inverted,
# so a relatively long window AND a short step make for best results

s = soundgen(samplingRate = samplingRate, addSilence = 50)
spec = spectrogram(s, samplingRate = samplingRate,
  wn = wn, windowLength = windowLength, step = step,
  zp = 0,  # otherwise it changes the window length and messes up istft
  padWithSilence = FALSE, output = 'original')

# Invert the spectrogram, attempting to guess the phase
# Note that we need to know the original windowLength, step, and wn
# (i.e., you have to know how the spectrogram was created)
s_new = invertSpectrogram(spec, samplingRate = samplingRate,
  windowLength = windowLength, step = step, wn = wn,
  initialPhase = 'spsi', nIter = 50, play = FALSE)

# Verify the quality of audio reconstruction
# playme(s, samplingRate); playme(s_new, samplingRate)

## Not run: 
# to improve the quality of reconstruction, increase the number of iterations
s_new = invertSpectrogram(spec, samplingRate = samplingRate,
  windowLength = windowLength, step = step, wn = wn,
  initialPhase = 'spsi', nIter = 500, play = FALSE)
playme(s, samplingRate); playme(s_new, samplingRate)
spectrogram(s, samplingRate)
spectrogram(s_new, samplingRate)

## End(Not run)

Match soundgen pars (experimental)

Description

Attempts to find settings for soundgen that will reproduce an existing sound. The principle is to mutate control parameters, trying to improve fit to target. The currently implemented optimization algorithm is simple hill climbing. Disclaimer: this function is experimental and may or may not work for particular tasks. It is intended as a supplement to - not replacement of - manual optimization. See https://cogsci.se/soundgen/sound_generation.html and https://cogsci.se/soundgen/matching/matching.html for more information.

Usage

matchPars(
  target,
  samplingRate = NULL,
  pars = NULL,
  specFun = "melspec",
  specFun_pars = list(),
  init = NULL,
  probMutation = 0.25,
  stepVariance = 0.1,
  maxIter = 50,
  minExpectedDelta = 0.001,
  compareSounds_pars = list(),
  verbose = TRUE,
  play = FALSE
)

Arguments

Value

A list containing the history of parameters tried and their final values (pars).

Examples

## Not run: 
target = soundgen(sylLen = 600, pitch = c(300, 200),
                  rolloff = -20, play = TRUE, plot = TRUE)
# we hope to reproduce this sound

# Match pars based on acoustic analysis alone, without any optimization.
# This *MAY* match temporal structure, pitch, and stationary formants
m1 = matchPars(target = target,
               samplingRate = 16000,
               maxIter = 0,  # no optimization, only acoustic analysis
               verbose = TRUE)
cand1 = do.call(soundgen, c(m1$pars, list(
  temperature = 0.001, play = TRUE, plot = TRUE)))

# Try to improve the match by optimizing rolloff
# (this may take a few minutes to run, and the results may vary)
m2 = matchPars(target = target,
               samplingRate = 16000,
               pars = 'rolloff',
               maxIter = 100,
               verbose = TRUE)
# rolloff should be moving from default (-12) to target (-20):
lapply(m2$history, function(x) x$pars$rolloff)
cand2 = do.call(soundgen, c(m2$pars, list(play = TRUE, plot = TRUE)))

## End(Not run)

Mix signals at specific SNR

Description

Adds two partly overlapping vectors, such as two waveforms or "signal" and "noise", at specific RMS or loudness levels. The location at which vector 2 ("noise") is pasted is defined by insertionPoint. Algorithm: both vectors are padded with zeros to match in length, rescaled, added, and re-normalized. All NA's are converted to 0. addVectors is the former name from soundgen 2.x, preserved for back-compatibility.

Usage

mixAtSNR(
  signal,
  noise,
  insertionPoint = 1L,
  normalize = TRUE,
  SNR = NULL,
  type = c("rms", "RMS", "peak", "loudness"),
  loudness_pars = list()
)

addVectors(
  signal,
  noise,
  insertionPoint = 1L,
  normalize = TRUE,
  SNR = NULL,
  type = c("rms", "RMS", "peak", "loudness"),
  loudness_pars = list()
)

Arguments

Value

A numeric vector.

Examples

v1 = 1:6
v2 = rep(100, 3)
addVectors(v1, v2, insertionPoint = 5, normalize = FALSE)
addVectors(v1, v2, insertionPoint = -4, normalize = FALSE)
addVectors(v1, rep(100, 15), insertionPoint = -4, normalize = FALSE)
# note the asymmetry: insertionPoint refers to the first arg
addVectors(v1, v2, insertionPoint = -4, normalize = FALSE)

v3 = rep(100, 15)
addVectors(v1, v3, insertionPoint = -4, normalize = FALSE)
addVectors(v2, v3, insertionPoint = 7, normalize = FALSE)
addVectors(1:6, 3:6, insertionPoint = 3, normalize = FALSE)

# add noise to a sound at a particular SNR
signal = rep(c(rep(0, 500), cospi(2 * 440 * (1:1000) / 16000)), 2)
noise = rnorm(length(signal))
osc(mixAtSNR(signal, noise, SNR = 6), samplingRate = 16000)
osc(mixAtSNR(signal, noise, SNR = 6, type = 'peak'), samplingRate = 16000)
osc(mixAtSNR(signal, noise, SNR = 6, type = 'loudness',
  loudness_pars = list(samplingRate = 16000, SPL = 80)), samplingRate = 16000)
osc(mixAtSNR(signal, noise, SNR = 12), samplingRate = 16000)
osc(mixAtSNR(signal, noise, SNR = 0), samplingRate = 16000)

# no rescaling (SNR = NULL)
osc(mixAtSNR(signal, noise, SNR = NULL), samplingRate = 16000)
osc(mixAtSNR(signal, noise/4, SNR = NULL), samplingRate = 16000)

Modulation spectrum

Description

Produces a modulation spectrum of waveform(s) or audio file(s). It begins with a spectrogram-like time-frequency representation and analyzes the modulation of the envelope in each frequency band. If specFun = 'audSpec', the sound is passed through a bank of bandpass filters with audSpectrogram. If specFun = 'STFT', we begin with an ordinary spectrogram produced with a Short-Time Fourier Transform. If msType = '2D', the modulation spectrum is a 2D Fourier transform of the spectrogram-like representation, with temporal modulation along the X axis and spectral modulation along the Y axis. A good visual analogy is decomposing the spectrogram into a sum of ripples of various frequencies and directions. If msType = '1D', the modulation spectrum is a matrix containing 1D Fourier transforms of each frequency band in the spectrogram, so the result again has modulation frequencies along the X axis, but the Y axis now shows the frequency of each analyzed band.

Usage

modulationSpectrum(
  x,
  samplingRate = NULL,
  scale = NULL,
  from = NULL,
  to = NULL,
  msType = c("2D", "1D"),
  specFun = "STFT",
  specFun_pars = list(),
  amRes = 5,
  maxDur = 5,
  specMethod = c("meanSpectrum", "spectrum"),
  logSpec = FALSE,
  logMPS = FALSE,
  power = 1,
  normalize = TRUE,
  roughRange = NULL,
  roughMean = 100,
  roughSD = 8,
  roughMinFreq = 1,
  amRange = c(10, 200),
  fluctRange = c(0.25, 30),
  fluctMean = 4,
  fluctSD = 12,
  returnMS = TRUE,
  returnComplex = FALSE,
  summaryFun = c("mean", "sd"),
  output = "all",
  averageMS = FALSE,
  reportEvery = NULL,
  cores = 1,
  plot = TRUE,
  savePlots = FALSE,
  embed = FALSE,
  logWarpX = NULL,
  logWarpY = NULL,
  quantiles = c(0.5, 0.8, 0.9),
  kernelSize = 5,
  kernelSD = 0.5,
  colorTheme = "bw",
  col = NULL,
  main = NULL,
  xlab = "Hz",
  ylab = NULL,
  xlim = NULL,
  ylim = NULL,
  width = 900,
  height = 500,
  units = "px",
  res = NA,
  ...
)

Arguments

Details

To calculate roughness and fluctuation depth, the symmetric modulation spectrum is folded by averaging its ±AM halves. Roughness is calculated as the proportion of the modulation spectrum within roughRange of temporal modulation frequencies or as a weighted version thereof. Fluctuation is calculated in the same manner as roughness, but in a lower frequency range of ~4 Hz (0.25 to 30 Hz). The frequency of amplitude modulation (amMsFreq, Hz) is calculated as the highest peak of the folded AM function within amRange, and its purity (amMsPurity, dB) as the ratio of this peak to the average of the folded AM function around the peak within amRange. For relatively short and steady sounds, set amRes = NULL and analyze the entire sound. For longer sounds and when roughness or AM vary over time, set amRes to get multiple measurements over time, but inspect the MS visually to make sure the relevant modulation frequencies are resolved. For multiple inputs, such as a list of waveforms or a 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 averageMS = TRUE.

The full STFT or auditory spectrogram is computed once and then sliced into temporal chunks. Roughness, fluctuation strength, amMsFreq, and amMsPurity are calculated from the per-fragment modulation spectra before Gaussian smoothing and before any log-warping used for plotting. Log-warps specified via logWarpX and logWarpY affect the plot only; the returned matrices are not log-warped.

Value

A list with three top-level elements: $summary, $detailed, and $aggregated.

$summary is a dataframe summarizing roughness, amMsFreq, and amMsPurity for each input (one row per file), or NULL if summaryFun = NULL.

$detailed contains detailed statistics per sound. If multiple sounds are analyzed, it is a named list of per-sound lists. If a single sound is analyzed, it is simplified to a single list. Components are selected via the output argument and may include:

original

averaged modulation spectrum across fragments, after logMPS, power, and normalize if requested, but before Gaussian smoothing and plotting transforms. Colnames are temporal modulation frequencies (Hz). Rownames are spectral modulation frequencies if msType = '2D' and frequencies of filters or spectrogram bands if msType = '1D'.

original_list

list of per-fragment modulation spectra.

modulation_spectrogram

a spectrogram-like representation showing how the frequency-averaged modulation spectrum changes over time; NA if there is only one fragment.

processed

modulation spectrum after Gaussian smoothing. This matrix is not log-warped; log-warping is applied only during plotting.

complex

complex modulation spectrum, returned only if returnComplex = TRUE, returnMS = TRUE, and msType = '2D'.

roughness

proportion of the modulation spectrum within roughRange, or a weighted version thereof if roughMean and roughSD are supplied, in percent. This is a vector if the sound is analyzed in multiple fragments, otherwise a single number.

roughness_list

a list containing frequencies, band-wise roughness values, and band amplitudes for each fragment.

roughness_spectrogram

a spectrogram-like matrix showing roughness per frequency band and per fragment.

fluctuation

proportion of the modulation spectrum within fluctRange, or a weighted version thereof if fluctMean and fluctSD are supplied, in percent. This is a vector if the sound is analyzed in multiple fragments, otherwise a single number.

fluctuation_list

a list containing frequencies, band-wise fluctuation values, and band amplitudes for each fragment.

fluctuation_spectrogram

a spectrogram-like matrix showing fluctuation strength per frequency band and per fragment.

amMsFreq

frequency of the highest AM peak within amRange. Like roughness, this can be a single number or a vector, depending on whether the sound is analyzed as a whole or in chunks.

amMsPurity

ratio of the AM peak at amMsFreq to the mean amplitude of the other folded AM values within amRange, in dB.

ampl

Root Mean Square amplitude of each analyzed fragment, divided by the internal audio scaling factor if present.

$aggregated is a list containing the original, processed, and (if requested) complex modulation spectra averaged across all successfully analyzed sounds. Only present if averageMS = TRUE and multiple sounds are analyzed.

References

• 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.

• Anikin, A. (2025) Acoustic estimation of voice roughness. Attention, Perception, & Psychophysics 87: 1771–1787.

See Also

plotMS spectrogram audSpectrogram analyze

Examples

s = soundgen(pitch = 440, amFreq = 100, amDep = 50)
ms = modulationSpectrum(s, samplingRate = 16000, amRes = NULL)
ms$detailed[c('roughness', 'amMsFreq', 'amMsPurity')]  # a single value for each
ms1 = modulationSpectrum(s, samplingRate = 16000, amRes = 5)
ms1$detailed[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, logMPS = TRUE,
  xlab = 'Temporal modulation, Hz', ylab = 'Spectral modulation, 1/kHz',
  colorTheme = 'matlab', main = 'Modulation spectrum', lty = 3)

# Plot a modulation spectrogram (the peak at 100 Hz shows AM)
spectrogram(s, 16000, specManual = ms$detailed$modulation_spectrogram,
  ylab = 'Modulation frequency, kHz', main = 'Modulation spectrogram')

# Plot a roughness spectrogram
spectrogram(s, 16000, specManual = ms$detailed$roughness_spectrogram,
  yScale = 'ERB', main = 'Roughness spectrogram')

# 1D instead of 2D
modulationSpectrum(s, 16000, msType = '1D', quantiles = NULL,
  col = soundgen:::jet.col(50), logWarpX = c(10, 2))

## Not run: 
# 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), plot = TRUE, yScale = 'ERB')
playme(s_long)
ms = modulationSpectrum(s_long, 16000)
# plot AM over time
plot(names(ms$detailed$amMsFreq), ms$detailed$amMsFreq,
  cex = ms$detailed$amMsPurity/5,
  xlab = 'Time, ms', ylab = 'AM frequency, Hz')
# plot roughness over time
spectrogram(s_long, 16000, ylim = c(0, 4),
  extraContour = list(x = ms$detailed$roughness / max(ms$detailed$roughness) * 4000, col = 'blue'))


# As with spectrograms, there is a tradeoff in time-frequency resolution
s = soundgen(pitch = 500, amFreq = 50, amDep = 100, sylLen = 500,
             samplingRate = 44100, plot = TRUE)
# playme(s, samplingRate = 44100)
ms = modulationSpectrum(s, samplingRate = 44100, amRes = NULL,
  specFun_pars = list(windowLength = 50, step = 50)) # poor temporal resolution
ms = modulationSpectrum(s, samplingRate = 44100, amRes = NULL,
  specFun_pars = list(windowLength = 5, step = 1))  # poor frequency resolution
ms = modulationSpectrum(s, samplingRate = 44100, amRes = NULL,
  specFun_pars = list(windowLength = 15, step = 3))  # a reasonable compromise

# Start with an auditory spectrogram instead of STFT
modulationSpectrum(s, 44100, specFun = 'audSpec', xlim = c(-100, 100))
modulationSpectrum(s, 44100, specFun = 'audSpec',
  logWarpX = c(10, 2), xlim = c(-500, 500),
  specFun_pars = list(nFilters = 32, filterType = 'gammatone', bandwidth = NULL))

# customize the plot
ms = modulationSpectrum(s, samplingRate = 44100,
  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
  logWarpX = c(10, 2),  # pseudo-log transform
  power = 2)                   # ^2
# Note the peaks at FM = 2/kHz (from "pitch = 500") and AM = 50 Hz (from
# "amFreq = 50")

# 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$detailed)  # separate MS and other descriptives per file
# ...or a single MS can be calculated by averaging across all files:
ms2 = modulationSpectrum('~/Downloads/temp', kernelSize = 11,
                         plot = FALSE, averageMS = TRUE)
plotMS(ms2$aggregated$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$detailed$roughness
# compare:
modulationSpectrum(s, samplingRate = 16000, maxDur = .5,
  xlim = c(-25, 25), colorTheme = 'seewave',
  power = 1)$detailed$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))$detailed$roughness
modulationSpectrum(s, samplingRate = 16000, logWarpX = NA)$detailed$roughness
modulationSpectrum(s, samplingRate = 16000, kernelSize = 17)$detailed$roughness

# Log-transform the spectrogram prior to 2D FFT (affects roughness):
modulationSpectrum(s, samplingRate = 16000, logSpec = FALSE)$detailed$roughness
modulationSpectrum(s, samplingRate = 16000, logSpec = TRUE)$detailed$roughness

# Use a lognormal weighting function to calculate roughness
# (instead of just % in roughRange)
modulationSpectrum(s, 16000, roughRange = NULL,
  roughMean = 75, roughSD = 3)$detailed$roughness
modulationSpectrum(s, 16000, roughRange = NULL,
  roughMean = 100, roughSD = 12)$detailed$roughness
# truncate weights outside roughRange
modulationSpectrum(s, 16000, roughRange = c(30, 150),
  roughMean = 100, roughSD = 1000)$detailed$roughness  # very large SD
modulationSpectrum(s, 16000, roughRange = c(30, 150),
  roughMean = NULL)$detailed$roughness  # same as above b/c SD --> Inf

# Complex modulation spectrum with phase preserved
ms = modulationSpectrum(soundgen(), samplingRate = 16000,
                        returnComplex = TRUE)
plotMS(log(abs(ms$detailed$complex + 1e-8)))  # note the symmetry
# compare:
plotMS(ms$detailed$original)

## End(Not run)

Morph sounds

Description

Takes two formulas for synthesizing two target sounds with soundgen and produces a number of intermediate forms (morphs), attempting to go from one target sound to the other in a specified number of equal steps. Normally you will want to set temperature very low; the tempEffects argument is not supported. Undefined values are set to reasonably intelligent defaults, but it is safer to specify the same parameters in both formulas. For instance, if formula1 includes "noise" and "mouth", set these parameters explicitly in formula2 as well.

Usage

morph(
  formula1,
  formula2,
  nMorphs,
  play = FALSE,
  savePath = NA,
  samplingRate = 16000
)

Arguments

Value

A list containing two lists, ($formulas and $sounds), each of length nMorphs. For ex., the formula for the second hybrid is m$formulas[[2]], and the waveform is m$sounds[[2]]

See Also

soundgen

Examples

## Not run: 
# write two formulas or copy-paste them from soundgen_app() or presets:
playback = c(TRUE, FALSE)[1]
# [a] to barking
m = morph(formula1 = list(repeatBout = 2),
          # equivalently: formula1 = 'soundgen(repeatBout = 2)',
          formula2 = presets$Misc$Dog_bark,
          nMorphs = 5, play = playback)
 # use $formulas to access formulas for each morph, $sounds for waveforms
 # m$formulas[[4]]
 # playme(m$sounds[[3]])

# morph intonation and vowel quality
m = morph(
  'soundgen(pitch = c(300, 250, 400),
            formants = c(350, 2900, 3600, 4700))',
  'soundgen(pitch = c(300, 700, 500, 300),
            formants = c(800, 1250, 3100, 4500))',
  nMorphs = 5, play = playback
)

# from a grunt of disgust to a moan of pleasure
m = morph(
  formula1 = 'soundgen(sylLen = 180, pitch = c(160, 160, 120), rolloff = -12,
    nonlinBalance = 70, subDep = 15, jitterDep = 2,
    formants = c(550, 1200, 2100, 4300, 4700, 6500, 7300),
    noise = data.frame(time = c(0, 180, 270), value = c(-20, -20, -30)),
    rolloffNoise = 0)',
  formula2 = 'soundgen(sylLen = 320, pitch = c(340, 330, 300),
    rolloff = c(-18, -16, -30), ampl = c(0, -10), formants = c(950, 1700, 3700),
    noise = data.frame(time = c(0, 300, 440), value = c(-25, -15, -35)),
    mouth = c(.4, .5), rolloffNoise = -5, attackLen = 30)',
  nMorphs = 8, play = playback
)

# from scream_010 to moan_515b
# (see online demos at http://cogsci.se/soundgen/humans/humans.html)
m = morph(
  formula1 = "soundgen(
    sylLen = 490,
    pitch = list(time = c(0, 80, 250, 370, 490),
    value = c(1000, 2900, 3200, 2900, 1000)),
    rolloff = c(-5, 0, -25), rolloffKHz = 0,
    temperature = 0.001,
    jitterDep = c(.5, 1, 0), shimmerDep = c(5, 15, 0),
    formants = c(1100, 2300, 3100, 4000, 5300, 6200),
    mouth = c(.3, .5, .6, .5, .3))",
  formula2 = "soundgen(sylLen = 520,
    pitch = c(300, 310, 300),
    ampl = c(0, -30),
    temperature = 0.001, rolloff = c(-18, -25),
    jitterDep = .05, shimmerDep = 2,
    formants = list(f1 = c(700, 900),
      f2 = c(1600, 1400),
      f3 = c(3600, 3500), f4 = c(4300, 4200)),
    mouth = c(.5, .3),
    noise = data.frame(time = c(0, 400, 660),
                       value = c(-40, -20, -70)),
    rolloffNoise = c(-5, -15))",
 nMorphs = 5, play = playback
)

## End(Not run)

Naive Bayes

Description

An implementation of a Naive Bayes classifier adapted to autocorrelated time series such as the type of nonlinear vocal phenomena in consecutive audio frames. All predictors must be continuous, and the outcome must be categorical. Cases with missing values are not deleted because the posterior probabilities of each outcome class can be calculated from different combinations of predictors on a case-by-case basis. Two optional modifications of a standard Naive Bayes algorithm can be made: (1) classifications can be "clumped" at the final stage, ensuring that every run or "epoch" of a particular predicted class is at least minLength steps long, and (2) priors can be continuously adapted based on static posterior probabilities of the preceding wlPrior observations if prior = 'dynamic'. naiveBayes_train calculates conditional means and standard deviations per class and global proportions of each class in the dataset. This is mostly useful because the output can be passed on to naiveBayes to save time if naiveBayes() is called in a loop with the same training dataset.

Usage

naiveBayes(
  formula,
  train,
  test = train,
  prior = c("static", "dynamic", "flat"),
  wlPrior = 3,
  wlClumper = NULL,
  runBack = TRUE,
  plot = FALSE
)

naiveBayes_train(formula, train)

Arguments

Value

A dataframe with new columns: "pr" = the predicted class membership, "[outcome][class]" = posterior probabilities per class, "like_[class]" = normalized likelihoods, "prior_[class]" = normalized priors, "priorF_[class]" / "priorB_[class]" = forward / backward normalized priors (only filled for prior = 'dynamic').

Functions

• naiveBayes_train(): Train a naive Bayes classifier

Examples

set.seed(151)
## create some fake data
df = data.frame(group = rep(c(
    rep('A', 150), rep('B', 50), rep('A', 120),
    rep('A', 100), rep('B', 30), rep('A', 90)
  ), 3))
df$group = as.factor(df$group)
df$x1 = rnorm(nrow(df), mean = ifelse(df$group == 'A', 3, 6), sd = 2)
df$x2 = rnorm(nrow(df), mean = ifelse(df$group == 'A', 2, -1), sd = 2)
boxplot(x1 ~ group, df)
boxplot(x2 ~ group, df)

## train the classifier
mod_train = naiveBayes_train(group ~ x1 + x2, train = df)
mod_train

## test on new data generated by the same process
test = data.frame(group = rep(c(
  rep('A', 90), rep('B', 40), rep('A', 150),
  rep('B', 40), rep('A', 130), rep('B', 30)
), 2))
test$group = as.factor(test$group)
test$x1 = rnorm(nrow(test), mean = ifelse(test$group == 'A', 3, 6), sd = 2)
test$x2 = rnorm(nrow(test), mean = ifelse(test$group == 'A', 2, -1), sd = 2)

# flat priors (same prior probability for each class)
nb_flat = naiveBayes(group ~ x1 + x2, train = mod_train, test = test,
  prior = 'flat', plot = TRUE)
# same as passing 'train' directly to the model, w/o calling naiveBayes_train():
nb_flat = naiveBayes(group ~ x1 + x2, train = df, test = test, prior = 'flat')
table(nb_flat$group, nb_flat$pr)
mean(nb_flat$group == nb_flat$pr)

# static priors (use original class proportions as prior class probabilities)
nb_static = naiveBayes(group ~ x1 + x2, train = mod_train, test = test,
  prior = 'static', wlClumper = NULL, plot = TRUE)
table(nb_static$group, nb_static$pr)
mean(nb_static$group == nb_static$pr)

# specify custom static priors
mod_train2 = mod_train
mod_train$table
mod_train2$table = list(A = .1, B = .9)  # sum to 1
nb_static2 = naiveBayes(group ~ x1 + x2, train = mod_train2, test = test,
  prior = 'static', wlClumper = NULL, plot = TRUE)
mean(nb_static2$group == nb_static2$pr)

# if we expect autocorrelation, ie class X is more likely a priori if the
# last few observations were also likely to be class X, we can use dynamic
# priors and/or clump the predicted classes (the latter imposes strong
# constraints on the predictions, but may be worth it if the data is known to
# be strongly "clumped", ie if we know classes occur in long'ish runs)
nb1 = naiveBayes(group ~ x1 + x2, train = mod_train, test = test,
  prior = 'dynamic', wlPrior = 10, plot = TRUE)
table(nb1$group, nb1$pr)
mean(nb1$group == nb1$pr)

nb2 = naiveBayes(group ~ x1 + x2, train = mod_train, test = test,
  prior = 'static', wlClumper = 10, plot = TRUE)
table(nb2$group, nb2$pr)
mean(nb2$group == nb2$pr)

nb3 = naiveBayes(group ~ x1 + x2, train = mod_train, test = test,
  prior = 'dynamic', wlPrior = 10, wlClumper = 10, plot = TRUE)
table(nb3$group, nb3$pr)
mean(nb3$group == nb3$pr)

# NAs in the data are not a problem
test1 = test
test1$x1[sample(1:nrow(test1), 100)] = NA
test1$x2[sample(1:nrow(test1), 10)] = NA
summary(test1)

nb4 = naiveBayes(group ~ x1 + x2, train = mod_train, test = test1,
  prior = 'dynamic', wlPrior = 10, plot = TRUE)
table(nb4$group, nb4$pr, useNA = 'ifany')
mean(nb4$group == nb4$pr, na.rm = TRUE)

Noise removal

Description

Removes noise by log-spectral subtraction. If a recording is affected by a steady noise with a relatively stable amplitude and spectrum (e.g., microphone hiss, crickets, MRI buzz, etc.), its log-spectrum can be simply subtracted from that of the signal. Algorithm: STFT to produce a log-spectrogram, subtract normalized noise spectrum, iSTFT to reconstitute the signal. Most of the work is done by addFormants.

Usage

noiseRemoval(
  x,
  samplingRate = NULL,
  from = NULL,
  to = NULL,
  noise,
  dB = 6,
  specificity = 1,
  windowLength = 50,
  step = NULL,
  overlap = 75,
  dynamicRange = 120,
  normalize = c("orig", "max", "none"),
  reportEvery = NULL,
  cores = 1,
  play = FALSE,
  saveAudio = FALSE,
  plot = FALSE,
  savePlots = FALSE,
  embed = FALSE,
  width = 900,
  height = 500,
  units = "px",
  res = NA,
  ...
)

Arguments

Value

The denoised audio as a numeric vector (multiple inputs return a list).

See Also

addFormants

Examples

s = soundgen(noise = list(time = c(-100, 400), value = -10),
  formantsNoise = list(f1 = list(freq = 3000, width = 25)),
  addSilence = 50, temperature = .001, plot = TRUE)
# Option 1: use part of the recording as noise profile
s1 = noiseRemoval(s, samplingRate = 16000, noise = c(0.05, 0.15),
  dB = 40, plot = TRUE)

## Not run: 
# Option 2: use a separate recording as noise profile
noise = soundgen(pitch = NA, noise = 0,
  formantsNoise = list(f1 = list(freq = 3000, width = 25)))
spectrogram(noise, 16000)
s2 = noiseRemoval(s, samplingRate = 16000, noise = noise,
  dB = 40, plot = TRUE)

# Option 3: provide noise spectrum as a matrix
spec_noise = spectrogram(
      noise, samplingRate = 16000,
      output = 'original', plot = FALSE)
s3 = noiseRemoval(s, samplingRate = 16000, noise = spec_noise,
  dB = 40, plot = TRUE)

# play with gain and specificity
s4 = noiseRemoval(s, samplingRate = 16000, noise = c(0.05, 0.15),
  dB = 60, specificity = 2, plot = TRUE)

# remove noise only from a section of the audio
s5 = noiseRemoval(s, samplingRate = 16000, from = .3, to = .5,
  noise = c(0.05, 0.15), dB = 60, plot = TRUE)

## End(Not run)

Normalize folder

Description

Normalizes the amplitude of all wav/mp3 files in a folder based on their peak or RMS amplitude or subjective loudness. This is good for playback experiments, which require that all sounds should have similar intensity or loudness. If preserveRelativeDif is set to TRUE, all recordings in the target folder get a boost in amplitude, so the loudest one reaches maxAmp, but the relative differences between the recordings are preserved (i.e., all files are boosted by the same amount - as much as possible to avoid clipping the loudest one).

Usage

normalizeFolder(
  myfolder,
  type = c("peak", "rms", "loudness"),
  maxAmp = 0,
  summaryFun = "mean",
  preserveRelativeDif = FALSE,
  windowLength = 50,
  step = NULL,
  overlap = 70,
  killDC = FALSE,
  windowDC = 200,
  cores = 1,
  savePath = NULL,
  reportEvery = NULL
)

Arguments

Details

Algorithm: first all files are rescaled to have the same peak amplitude of maxAmp dB. If type = 'peak', the process ends here. If type = 'rms', there are two additional steps. First the original RMS amplitude of all files is calculated per frame by getRMS. The "quietest" sound with the lowest summary RMS value is not modified, so its peak amplitude remains maxAmp dB. All the remaining sounds are rescaled linearly, so that their summary RMS values becomes the same as that of the "quietest" sound, and their peak amplitudes become smaller, <maxAmp. Finally, if type = 'loudness', the subjective loudness of each sound is estimated by getLoudness, which assumes frequency sensitivity typical of human hearing. The following normalization procedure is similar to that for type = 'rms'. NB: because loudness is not a simple linear function of SPL, loudness normalization is only approximate; reiterate the process several times to improve the precision of loudness normalization.

Value

Does not return anything, only saves the normalized audio files.

See Also

getRMS analyze getLoudness

Examples

## Not run: 
# put a few short audio files in a folder, eg '~/Downloads/temp'
target = '~/Downloads/temp'
save_in_folder = paste0(target, '/normalized')
getRMS(target, summaryFun = 'mean')$summary  # different
normalizeFolder(target, type = 'rms', summaryFun = 'mean',
  savePath = save_in_folder)
getRMS(save_in_folder, summaryFun = 'mean')$summary  # same
# If the saved audio files are treated as stereo with one channel missing,
# try reconverting with ffmpeg (saving is handled by tuneR::writeWave)

## End(Not run)

Oscillogram

Description

Plots the oscillogram (waveform) of a sound on a linear or logarithmic scale (in dB). To get a dB scale, centers and normalizes the sound, then takes a logarithm of the positive part and a flipped negative part, which is analogous to "Waveform (dB)" view in Audacity. For more plotting options, check the "oscillo" function in the seewave package.

Usage

osc(
  x,
  samplingRate = NULL,
  scale = NULL,
  from = NULL,
  to = NULL,
  dynamicRange = 80,
  dB = FALSE,
  returnWave = FALSE,
  reportEvery = NULL,
  cores = 1,
  plot = TRUE,
  savePlots = FALSE,
  embed = FALSE,
  main = NULL,
  xlab = NULL,
  ylab = NULL,
  ylim = NULL,
  bty = "n",
  midline = TRUE,
  maxPoints = 10000,
  width = 900,
  height = 500,
  units = "px",
  res = NA,
  ...
)

Arguments

Value

If returnWave = TRUE, returns the waveform as a numeric vector on the original or dB scale; otherwise, does not return anything.

Examples

sound = sin(1:2000/10) * approx(c(1, .01, .5), n = 2000)$y

# Oscillogram on a linear scale without bells and whistles, just base R
plot(sound, type = 'l')

# Oscillogram options with soundgen
osc(sound)             # linear
osc(sound, dB = TRUE)  # dB

# For numeric vectors, indicate samplingRate and scale (max amplitude)
osc(sound, samplingRate = 1000, scale = 100, dB = TRUE)

# Embellish and customize the plot
o = osc(sound, samplingRate = 1000, dB = TRUE, midline = FALSE,
        main = 'My waveform', col = 'blue', returnWave = TRUE)
abline(h = -80, col = 'orange', lty = 3)
o[1:10]  # the waveform in dB

# Wave object
data(speechEx, package = 'soundgen')
osc(speechEx, from = 2, to = 2.5, dB = TRUE)

# for long files, reduce the resolution to plot quickly (careful: if the
# resolution is too low, antialiasing may cause artifacts)
osc(speechEx, dB = TRUE, maxPoints = 2500)

## Not run: 
# saves oscillograms of all audio files in a folder
osc('~/Downloads/temp', savePlots = TRUE)

## End(Not run)

Phasegram

Description

Produces a phasegram of a sound or another time series, which is a collection of Poincare sections cut through phase portraits of consecutive frames. The x axis is time, just as in a spectrogram, the y axis is a slice through the phase portrait, and the color shows the density of trajectories at each point of the phase portrait.

Usage

phasegram(
  x,
  samplingRate = NULL,
  from = NULL,
  to = NULL,
  windowLength = 10,
  step = NULL,
  overlap = 50,
  timeLag = NULL,
  theilerWindow = NULL,
  nonlinStats = c("ed", "d2", "ml", "sur"),
  ed_pars = list(max.embedding.dim = 15),
  d2_pars = list(min.embedding.dim = 2, min.radius = 0.001, n.points.radius = 20),
  ml_pars = list(min.embedding.dim = 2, radius = 0.001),
  sur_pars = list(FUN = nonlinearTseries::timeAsymmetry, K = 20),
  bw = 0.01,
  bins = 5/bw,
  reportEvery = NULL,
  cores = 1,
  rasterize = FALSE,
  plot = TRUE,
  savePlots = FALSE,
  embed = FALSE,
  colorTheme = "bw",
  col = NULL,
  xlab = "Time",
  ylab = "",
  main = NULL,
  width = 900,
  height = 500,
  units = "px",
  res = NA,
  ...
)

Arguments

Details

Algorithm: the input sound is normalized to [-1, 1] and divided into consecutive frames windowLength ms long without multiplying by any windowing function (unlike in STFT). For each frame, a phase portrait is obtained by time-shifting the frame by timeLag ms. A Poincare section is taken through the phase portrait (currently at a fixed angle, namely the default in poincareMap), giving the intersection points of trajectories with this bisecting line. The density of intersections is estimated with a smoothing kernel of bandwidth bw (as an alternative to using histogram bins). The density distributions per frame are stacked together into a phasegram (output: orig). The density values in orig are normalized by the global maximum across all frames. The resulting phasegram can optionally be rasterized to smooth it for plotting (output: rasterized); the rasterized matrix is additionally normalized row by row for display.

Value

For a single input, a list of three components:

orig

the full phasegram as a data frame. $time is the middle of each frame (ms), $x is the coordinate along the Poincare section (approximately on the normalized audio scale), and $y is the density of intersections of system trajectories with the Poincare section, normalized by the global maximum across all frames. Failed or flat frames are represented by NA rows

rasterized

the rasterized phasegram as a numeric matrix, or NULL if rasterize = FALSE. Rows correspond to time frames and columns correspond to bins along the Poincare-section coordinate. Values are row-normalized for display, so each non-empty row has a maximum of 1. If no valid Poincare intersections are found, a zero-valued matrix is returned

descriptives

per-frame descriptives as a data frame. Always included are time (ms), shannon = normalized Shannon entropy of Poincare sections, and nPeaks = log-normalized number of peaks in the density distribution of Poincare sections. If requested via nonlinStats, also includes ed = optimal number of embedding dimensions, d2 = correlation dimension, ml = maximum Lyapunov exponent (positive values suggest chaos), and sur = stochasticity index from surrogate data testing, rescaled to approximately 0 = deterministic and 1 = stochastic

For multiple inputs, a list of such per-input results is returned.

References

• Herbst, C. T., Herzel, H., Švec, J. G., Wyman, M. T., & Fitch, W. T. (2013). Visualization of system dynamics using phasegrams. Journal of the Royal Society Interface, 10(85), 20130288.

• Huffaker, R., Huffaker, R. G., Bittelli, M., & Rosa, R. (2017). Nonlinear time series analysis with R. Oxford University Press.

Examples

target = soundgen(sylLen = 300, pitch = c(350, 420, 420, 410, 340) * 3,
  subDep = c(0, 0, 60, 50, 0, 0) / 2, addSilence = 0, plot = TRUE)
# Nonlinear statistics are also returned (slow - disable by setting
# nonlinStats = NULL if these are not needed)
ph = phasegram(target, 16000, nonlinStats = NULL)

## Not run: 
ph = phasegram(target, 16000, windowLength = 20, step = 20,
  rasterize = TRUE, bw = .01, bins = 150)
ph$descriptives

# Unfortunately, phasegrams are greatly affected by noise. Compare:
target2 = soundgen(sylLen = 300, pitch = c(350, 420, 420, 410, 340) * 3,
  subDep = c(0, 0, 60, 50, 0, 0), noise = -30, jitterDep = .4,
  rolloff = -5, addSilence = 0, plot = TRUE)
ph2 = phasegram(target2, 16000, nonlinStats = NULL)

# low-pass filtering may help a bit
target2_lowpass = bandpass(target2, 16000, upr = 2500)
phasegram(target2_lowpass, 16000, nonlinStats = NULL)

s2 = soundgen(sylLen = 3000, addSilence = 0, temperature = 1e-6,
  pitch = c(380, 550, 500, 220), subDep = c(0, 0, 40, 0, 0, 0, 0, 0),
  amDep = c(0, 0, 0, 0, 80, 0, 0, 0), amFreq = 80,
  jitterDep = c(0, 0, 0, 0, 0, 3), plot = TRUE, yScale = 'bark')
phasegram(s2, 16000, windowLength = 10, nonlinStats = NULL, bw = .001)
phasegram(s2, 16000, windowLength = 10, nonlinStats = NULL, bw = .02)

## End(Not run)

Pitch descriptives

Description

Provides common descriptives of time series such as pitch contours, including measures of average / range / variability / slope / inflections etc. Several degrees of smoothing can be applied consecutively. The summaries are produced on the original and log-transformed scales, so this is meant to be used on frequency-related variables in Hz.

Usage

pitchDescriptives(
  x,
  step = NULL,
  timeUnit = NULL,
  smoothBW = c(NA, 10, 1),
  inflThres = 0.2,
  summaryFun = c("mean", "sd"),
  extraSummaryFun = c(),
  ref = 16.3516,
  ptvStep = NULL,
  ptvTime = 0.5,
  ptvFreq = 1,
  plot = FALSE
)

Arguments

Value

A list with three elements: summary (a dataframe with columns containing summaries of one or multiple inputs, one input per row), syllables (a dataframe or list of dataframes with syllable timings, where a syllable is a contiguous non-NA run of the pitch contour), and ptv (a dataframe or list of dataframes with the instantaneous PTV contour). The descriptives in summary are as follows:

duration

total duration, s

durDefined

duration after omitting leading and trailing NAs

propDefined

percentage of input with non-NA value, eg percentage of voiced frames if the input is pitch

start, start_oct, end, end_oct

the first and last values on the original scale and in octaves above C0 (16.3516 Hz)

mean, median, max, min

average and extreme values on the original scale

mean_oct, median_oct, min_oct, max_oct

same in octaves above C0

time_max, time_min

the location of minimum and maximum relative to durDefined, 0 to 1

range, range_sem, sd, sd_sem

range and standard deviation on the original scale and in semitones

CV

coefficient of variation = sd/mean (provided for historical reasons)

meanSlope, meanSlope_sem

mean slope in Hz/s or semitones/s (NB: does not depend on duration or missing values)

meanAbsSlope, meanAbsSlope_sem

mean absolute slope (modulus, ie rising and falling sections no longer cancel out)

maxAbsSlope, maxAbsSlope_sem

the steepest slope

Examples

x = c(NA, NA, 405, 441, 459, 459, 460, 462, 462, 458, 458, 445, 458, 451,
444, 444, 430, 416, 409, 403, 403, 389, 375, NA, NA, NA, NA, NA, NA, NA, NA,
NA, 183, 677, 677, 846, 883, 886, 924, 938, 883, 946, 846, 911, 826, 826,
788, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 307,
307, 368, 377, 383, 383, 383, 380, 377, 377, 377, 374, 374, 375, 375, 375,
375, 368, 371, 374, 375, 361, 375, 389, 375, 375, 375, 375, 375, 314, 169,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 238, 285, 361, 374, 375, 375,
375, 375, 375, 389, 403, 389, 389, 375, 375, 389, 375, 348, 361, 375, 348,
348, 361, 348, 342, 361, 361, 361, 365, 365, 361, 966, 966, 966, 959, 959,
946, 1021, 1021, 1026, 1086, 1131, 1131, 1146, 1130, 1172, 1240, 1172, 1117,
1103, 1026, 1026, 966, 919, 946, 882, 832, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA)
plot(x, type = 'b')
ci95 = function(x) diff(quantile(na.omit(x), probs = c(.025, .975)))
pd = pitchDescriptives(
  x, step = .025,
  smoothBW = c(NA, 10, 1),   # original + smoothed at 10 Hz and 1 Hz
  inflThres = c(NA, .2, .2), # different for each level of smoothing
  extraSummaryFun = 'ci95',  # user-defined, here 95% coverage interval
  plot = TRUE
)
pd

## Not run: 
# a single file
data(speechEx, package = 'soundgen')
a = analyze(speechEx)
pd1 = pitchDescriptives(a$detailed[, c('time', 'pitch')],
                        inflThres = NA, plot = TRUE)
pd2 = pitchDescriptives(a$detailed[, c('time', 'pitch')],
                        inflThres = c(0.1, 0.1, .5), plot = TRUE)

# multiple files returned by analyze()
an = analyze('~/Downloads/temp')
pd = pitchDescriptives(an$detailed)
pd

## End(Not run)

Pitch smoothing as in Praat

Description

Smooths an intonation (pitch) contour with a low-pass filter, as in Praat (http://www.fon.hum.uva.nl/praat/). Algorithm: interpolates missing values (voiceless frames), performs FFT to obtain the spectrum, multiplies by a Gaussian filter, performs an inverse FFT, and fills the missing values back in. The bandwidth parameter controls the width of the low‑pass filter; larger values retain more high‑frequency detail. Similar to a special case of low-pass filtering with bandpass, but faster and can handle NAs.

Usage

pitchSmoothPraat(
  pitch,
  bandwidth,
  samplingRate,
  plot = FALSE,
  preprocess = TRUE
)

Arguments

Value

Numeric vector of the same length as input.

See Also

analyze

Examples

pitch = c(NA, NA, 405, 441, 459, 459, 460, 462, 462, 458, 458, 445, 458, 451,
444, 444, 430, 416, 409, 403, 403, 389, 375, NA, NA, NA, NA, NA, NA, NA, NA,
NA, 183, 677, 677, 846, 883, 886, 924, 938, 883, 946, 846, 911, 826, 826,
788, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 307,
307, 368, 377, 383, 383, 383, 380, 377, 377, 377, 374, 374, 375, 375, 375,
375, 368, 371, 374, 375, 361, 375, 389, 375, 375, 375, 375, 375, 314, 169,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 238, 285, 361, 374, 375, 375,
375, 375, 375, 389, 403, 389, 389, 375, 375, 389, 375, 348, 361, 375, 348,
348, 361, 348, 342, 361, 361, 361, 365, 365, 361, 966, 966, 966, 959, 959,
946, 1021, 1021, 1026, 1086, 1131, 1131, 1146, 1130, 1172, 1240, 1172, 1117,
1103, 1026, 1026, 966, 919, 946, 882, 832, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA)
pitchSmoothPraat(pitch, bandwidth = 10, samplingRate = 40, plot = TRUE)
pitchSmoothPraat(pitch, bandwidth = 2, samplingRate = 40, plot = TRUE)

Interactive pitch tracker

Description

Starts a shiny app for manually editing pitch contours. The settings in the panels on the left correspond to arguments to analyze - see analyze and https://cogsci.se/soundgen/acoustic_analysis.html for help and examples. You can verify the pitch contours first, and then feed them back into analyze (see examples). Use hotkeys (eg spacebar to play/stop) and avoid working with very large files.

Usage

pitch_app(...)

Arguments

Value

A list with the last used settings ($settings) plus the output of analyze() for each file from the last file queue with two additional columns: "time" and "pitch". NB: only the results of the most recent file queue are returned, so don't press "Load audio" repeatedly if you need the output returned to R (the csv file with results should still be saved correctly). When proceeding to the next file in the queue, the app saves to disk a backup .csv file with one row per audio file. When the orange "Download results" button is clicked, a context menu pops up offering to terminate the app - if that happens, the results are also returned directly into R. To process pitch contours further in R, work directly with my_pitch[[myfile]]$detailed$time and my_pitch[[myfile]]$detailed$pitch or, if loading the csv file, do something like:

a = read.csv('~/Downloads/output.csv', stringsAsFactors = FALSE)
pitch = as.numeric(unlist(strsplit(a$pitch, ',')))
mean(pitch, na.rm = TRUE); sd(pitch, na.rm = TRUE)

Suggested workflow

Start by setting the basic analysis settings such as pitchFloor, pitchCeiling, silence, etc. Then click "Load audio" to upload one or several audio files (wav/mp3). Long files will be very slow, so please cut your audio into manageable chunks (ideally <10 s). If Shiny complains that maximum upload size is exceeded, you can increase it, say to 30 MB, with options(shiny.maxRequestSize = 30 * 1024^2). Once the audio has been uploaded to the browser, fine-tune the analysis settings as needed, edit the pitch contour in the first file to your satisfaction, then click "Next" to proceed to the next file, etc. Remember that setting a reasonable prior is often faster than adjusting the contour one anchor at a time. When done, click "Save results". If working with many files, you might want to save the results occasionally in case the app crashes (although you should still be able to recover your data if it does - see below).

How to edit pitch contours

Left-click to add a new anchor, double-click to remove it or unvoice the frame. Each time you make a change, the entire pitch contour is re-fit, so making a change in one frame can affect the path through candidates in adjacent frames. You can control this behavior by changing the settings in Out/Path and Out/Smoothing. If correctly configured, the app corrects the contour with only a few manual values - you shouldn't need to manually edit every single frame. For longer files, you can zoom in/out and navigate within the file. You can also select a region to voice/unvoice or shift it as a whole or to set a prior based on selected frequency range.

Recovering lost data

Every time you click "next" or "last" to move in between files in the queue, the output you've got so far is saved in a temporary backup file. If the app crashes or is closed without saving the results, this backup file preserves your data. To recover it, restart pitch_app() - a dialog box will pop up and ask whether you want to append the old data to the new one. Even so, save your data regularly to be on the safe side!

See Also

formant_app annotation_app

Examples

## Not run: 
# Recommended workflow for analyzing a lot of short audio files
path_to_audio = '~/Downloads/temp'  # our audio lives here

# STEP 1: extract manually corrected pitch contours
my_pitch = pitch_app()  # runs in default browser such as Firefox or Chrome
# To change system default browser, run something like:
options('browser' = '/usr/bin/firefox')  # path to the executable on Linux

# You can pass presets with your preferred parameter values:
my_pitch = pitch_app(windowLength = 20, step = 10,
  pitchMethods = c('dom', 'autocor', 'cep'), spec_ylim = c(0, 6))

# Object "my_pitch" contains the output, notably the time-pitch matrix
plot(my_pitch[[1]]$detailed$time, my_pitch[[1]]$detailed$pitch, type = 'b',
  xlab = 'Time, ms', ylab = 'Pitch, Hz')

# Run the app with previously used settings
my_pitch2 = do.call(pitch_app, my_pitch$settings)

# save the complete output, including the settings used
saveRDS(my_pitch2, 'my_pitch_analysis.rds')

# STEP 2: run analyze() with manually corrected pitch contours to obtain
# accurate descriptives like the proportion of energy in harmonics above f0,
# etc. This also gives you formants and loudness estimates (disabled in
# pitch_app to speed things up)
df2 = analyze(
  path_to_audio,
  pitchMethods = 'autocor',  # needed for calculating HNR
  nFormants = 5,        # now we can measure formants as well
  pitchManual = my_pitch
  # or, if loading the output of pitch_app() from the disk:
  # pitchManual = '~/Downloads/output.csv'
  # pitchManual = '~/path_to_some_folder/my_pitch_contours.rds
)

# STEP 3: add other acoustic descriptors, for ex.
df3 = segment(path_to_audio)

# STEP 4: merge df2, df3, df4, ... in R or a spreadsheet editor to have all
# acoustic descriptives together

# To verify your pitch contours and/or edit them later, copy output.csv to
# the folder with your audio, run pitch_app(), and load the audio + csv
# together. The saved pitch contours are treated as manual anchors

## End(Not run)

Play audio

Description

Plays one or more sounds: wav/mp3 file(s), Wave objects, or numeric vectors. Recommended players on Linux: "play" from the "SoX" library (default), "aplay" (ALSA).

Usage

playme(x, samplingRate = NULL, player = NULL, from = NULL, to = NULL)

Arguments

Value

Does not return anything, only plays the audio.

Examples

## Not run: 
# Play an audio file:
data('speechEx', package = 'soundgen')
playme(speechEx, from = 1, to = 2)
playme('pathToMyAudio/audio.wav')

# Create and play a numeric vector:
f0_Hz = 440
sound = sin(2 * pi * f0_Hz * (1:16000) / 16000)
playme(sound, 16000)
playme(sound, 16000, from = .1, to = .5)  # play from 100 to 500 ms

# Set default player for the current R session
# "play" (SoX library, default on Linux)
options(soundgen.player = 'play')
# "afplay" (ALSA, default on MacOS)
options(soundgen.player = 'afplay')
# VLC path (default on Windows)
options(soundgen.player = "C:\\Program Files\\VideoLAN\\VLC\\vlc.exe")
# view currently set player
options()$soundgen.player

# Or specify player per call
playme(sound, player = "C:\\Program Files\\VideoLAN\\VLC\\vlc.exe")

## End(Not run)

Plot modulation spectrum

Description

Plots a single modulation spectrum returned by modulationSpectrum. The result is the same as the plot produced by modulationSpectrum, but calling plotMS is handy for processed modulation spectra - for instance, for plotting the difference between the modulation spectra of two sounds or groups of sounds.

Usage

plotMS(
  ms,
  X = as.numeric(colnames(ms)),
  Y = as.numeric(rownames(ms)),
  quantiles = c(0.5, 0.8, 0.9),
  colorTheme = "bw",
  col = NULL,
  logWarpX = NULL,
  logWarpY = NULL,
  main = NULL,
  xlab = "Hz",
  ylab = "1/kHz",
  xlim = NULL,
  ylim = NULL,
  audio = NULL,
  extraY = TRUE,
  ...
)

Arguments

Value

Does not return anything, just plots the MS.

Examples

ms1 = modulationSpectrum(runif(4000), samplingRate = 16000, plot = TRUE)
plotMS(ms1$detailed$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$detailed$original,
  nr = nrow(ms1$detailed$original), nc = ncol(ms1$detailed$original))
# plot the difference
plotMS(log(ms1$detailed$original / ms2_resized), quantiles = NULL,
  col = colorRampPalette(c('blue', 'yellow')) (50))

Presets

Description

A library of presets for easy generation of a few representative sounds with soundgen.

Usage

presets

Format

A list of length 4.

Prosody

Description

Exaggerates or flattens the intonation by performing a dynamic pitch shift, changing pitch excursion from its original median value without changing the formants. This is a particular case of pitch shifting, which is performed with shiftPitch. The result is likely to be improved if manually corrected pitch contours are provided. Depending on the nature of audio, the settings that control pitch shifting may also need to be fine-tuned with the shiftPitch_pars argument. Any NAs in pitch contour are treated as voiceless fragments, and their pitch is not modified. If the NAs are not actual silences / voiceless frames, they should be interpolated before passing the pitch contour to prosody().

Usage

prosody(
  x,
  samplingRate = NULL,
  multProsody,
  analyze_pars = list(),
  shiftPitch_pars = list(),
  pitchManual = NULL,
  reportEvery = NULL,
  cores = 1,
  play = FALSE,
  saveAudio = FALSE,
  plot = FALSE,
  savePlots = FALSE,
  embed = FALSE,
  width = 900,
  height = 500,
  units = "px",
  res = NA
)

Arguments

Value

If the input is a single audio (file, Wave, or numeric vector), returns the processed waveform as a numeric vector with the original sampling rate and scale. If the input is a folder with several audio files, returns a list of processed waveforms, one for each file.

See Also

shiftPitch

Examples

s = soundgen(sylLen = 200, pitch = c(150, 220), addSilence = 50,
             plot = TRUE, yScale = 'log')
# playme(s)
s1 = prosody(s, 16000, multProsody = 2,
  analyze_pars = list(windowLength = 30, step = 15),
  shiftPitch_pars = list(windowLength = 20, step = 5, freqWindow = 300),
  plot = TRUE)
# playme(s1)
# spectrogram(s1, 16000, yScale = 'log')

## Not run: 
data('speechEx', package = 'soundgen')
samplingRate = speechEx@samp.rate
spectrogram(speechEx, yScale = 'log', ylim = c(.05, 4))
# playme(speechEx)

# start with exaggerated prosody, then flat towards the end
speech1 = prosody(speechEx,
  multProsody = list(time = c(0, 1), value = c(1.5, 0.1)),
  analyze_pars = list(windowLength = 40, step = 10,
  pitchMethods = c('dom', 'autocor', 'cep', 'spec')),
  shiftPitch_pars = list(freqWindow = 400))
spectrogram(speech1, samplingRate, yScale = 'log')
# playme(speech1, samplingRate)

# process all audio files in a folder
s4 = prosody('~/Downloads/temp', multProsody = 2,
  savePlots = TRUE, saveAudio = TRUE)
str(s4)  # returns a list with audio (+ saves it to disk)

## End(Not run)

Report CI

Description

A simple function for formatting and printing estimates with their confidence intervals or similar statistics that follow a "X [Y, Z]" format. Takes a numeric vector or matrix with three elements / columns: estimate, lower quantile from a CI, and upper quantile from a CI. For each row, it prints the result as "estimate [CI]".

Usage

reportCI(n, digits = 2, suffix = NULL)

Arguments

Value

Does not return anything, just prints the estimate and CI.

Examples

n = rnorm(100)
reportCI(quantile(n, probs = c(.5, .025, .975)))

a = data.frame(fit = c(3, 5, 7),
               lwr = c(1, 4, 6.5),
               upr = c(5, 6, 7.1))
reportCI(a, 1)
reportCI(a, 1, ' cm')
reportCI(a, 1, '%, 95% CI')

Report time

Description

Provides a nicely formatted "estimated time left" in loops plus a summary upon completion.

Usage

reportTime(
  i,
  time_start,
  nIter = NULL,
  reportEvery = NULL,
  jobs = NULL,
  prefix = ""
)

Arguments

Value

Does not return anything, just prints the time left.

Examples

time_start = proc.time()
nIter = 100
for (i in 1:nIter) {
  Sys.sleep(i ^ 1.02 / 10000)
  reportTime(i, time_start, nIter,
    jobs = (1:100) ^ 1.02, prefix = 'Chain 1: ')
}

# Unknown number of iterations:
time_start = proc.time()
for (i in 1:20) {
  Sys.sleep(i ^ 2 / 10000)
  reportTime(i = i, time_start = time_start,
  jobs = (1:20) ^ 2, reportEvery = 5)
}

## Not run: 
# when analyzing a bunch of audio files, their size is a good estimate
# of how long each will take to process
time_start = proc.time()
filenames = list.files('~/Downloads/temp', pattern = "*.wav|.mp3",
  full.names = TRUE)
filesizes = file.info(filenames)$size
for (i in seq_along(filenames)) {
  # ...do what you have to do with each file...
  reportTime(i = i, time_start = time_start, nIter = length(filenames),
             jobs = filesizes)
}

## End(Not run)

Resample a vector

Description

Changes the sampling rate of a sound or another time series without aliasing. Similar in spirit and speed to resample in the signal package, but also works for short vectors and preserves NAs approximately after resampling. Use resample() for interpolation when the points are evenly spaced, when there are missing values, or when we care about aliasing. When the points are unevenly spaced, use interpolate instead, but then there should not be any missing values.

Usage

resample(
  x,
  samplingRate = NULL,
  samplingRate_new = NULL,
  mult = NULL,
  len = NULL,
  lowPass = TRUE,
  interpol = "splineFC",
  na.rm = FALSE,
  reportEvery = NULL,
  cores = 1,
  saveAudio = FALSE,
  plot = FALSE,
  savePlots = FALSE,
  embed = FALSE,
  width = 900,
  height = 500,
  units = "px",
  res = NA
)

Arguments

Details

Algorithm: to downsample, applies a low-pass filter, then decimates with spline; to upsample, performs a specified type of interpolation (defaults to monotone Fritsch‑Carlson splines, see interpolate), then applies a low-pass filter. NAs can be interpolated or preserved in the output. The length of output is determined, in order of precedence, by len / mult / samplingRate_new. For simple vector operations, this is very similar to approx, but the leading and trailing NAs are also preserved when na.rm = FALSE (see examples).

Value

A numeric vector.

See Also

interpolate interpolateNA

Examples

## Example 1: a short vector with NAs
x = c(NA, 1, 2, 3, NA, NA, 6, 9, 8, NA)

# downsample
print(resample(x, len = 5))  # NAs are propagated
print(resample(x, len = 5, na.rm = TRUE))  # NAs are interpolated
print(resample(x, mult = 1/2, na.rm = TRUE, plot = TRUE)) # same

# upsample
resample(x, mult = 3.5, lowPass = FALSE, plot = TRUE)  # just approx
resample(x, mult = 3.5, lowPass = TRUE, plot = TRUE) # low-pass + approx
resample(x, mult = 3.5, lowPass = FALSE, na.rm = TRUE, plot = TRUE)

# change the method of interpolation
resample(x, mult = 15, lowPass = FALSE, interpol = 'pchip', plot = TRUE)
resample(x, mult = 15, lowPass = FALSE,
  interpol = interpol_loess(span = .6), plot = TRUE)

## Example 2: a sound
silence = rep(0, 10)
samplingRate = 1000
fr = seq(100, 300, length.out = 400)
x = c(silence, sin(cumsum(fr) * 2 * pi / samplingRate), silence)
spectrogram(x, samplingRate)

# downsample
x1 = resample(x, mult = 1 / 2.5)
spectrogram(x1, samplingRate / 2.5)  # no aliasing
# cf:
x1bad = resample(x, mult = 1 / 2.5, lowPass = FALSE)
spectrogram(x1bad, samplingRate / 2.5)  # aliasing

# upsample
x2 = resample(x, mult = 3)
spectrogram(x2, samplingRate * 3)  # nothing above the old Nyquist
# cf:
x2bad = resample(x, mult = 3, lowPass = FALSE)
spectrogram(x2bad, samplingRate * 3)  # high-frequency artifacts

## Not run: 
# Example 3: resample all audio files in a folder to 8000 Hz
resample('~/Downloads/temp', saveAudio = TRUE,
         samplingRate_new = 8000, savePlots = TRUE)

## End(Not run)

Reverb & echo

Description

Adds reverberation and/or echo to a sound parametrically or based on the impulse response recorded in a particular real environment. Algorithm for parametric reverb: adds time-shifted copies of the signal weighted by a decay function. In simple terms: we specify how much and when the sound rebounds back (e.g., from the walls of a room) and add these time-shifted copies to the original - optionally, with some spectral filtering.

Usage

reverb(
  x,
  samplingRate = NULL,
  echoDelay = 200,
  echoLevel = -20,
  reverbDelay = 70,
  reverbSpread = 130,
  reverbLevel = -25,
  reverbDensity = 50,
  reverbType = "gaussian",
  IR = NULL,
  filter = list(),
  dynamicRange = 80,
  output = c("audio", "detailed"),
  play = FALSE,
  saveAudio = FALSE,
  reportEvery = NULL,
  cores = 1,
  seed = NULL
)

Arguments

Value

If output = "audio", returns the processed waveform as a numeric vector. If output = "detailed", returns a list: $rvb_win = Gaussian window used to create the reverberation, $rvb = reverberation effect, $echo = echo effect, $effect = both reverb and echo mixed together, $audio = processed audio.

Examples

data('speechEx', package = 'soundgen')
sr = speechEx@samp.rate
s = speechEx@left[1:16000]
# playme(s)

# parametric reverb
s_rev = reverb(s, sr)
# playme(s_rev)

# reverb of a particular real place based on its recorded IR
data('IR', package = 'soundgen')
s_ir = reverb(s, sr, IR = IR)
# playme(s_ir)

## Not run: 
# double echo, no reverb
s1 = reverb(s, sr, reverbLevel = NULL,
            echoDelay = c(250, 800), echoLevel = c(-15, -25))
playme(s1)
spectrogram(s1, 16000, ylim = c(0, 4))

# only reverb (indoors)
s2 = reverb(s, sr, echoDelay = NULL,
            reverbDelay = 70, reverbSpread = 130,
            reverbLevel = -20, reverbDensity = 20)
playme(s2)
spectrogram(s2, 16000, ylim = c(0, 4))

# reverb (caves)
s3 = reverb(s, sr, echoDelay = NULL,
            reverbDelay = 600, reverbSpread = 1500,
            reverbLevel = -10, reverbDensity = 100)
playme(s3)
spectrogram(s3, 16000, ylim = c(0, 4))

# both echo and reverb with high frequencies emphasized
s4 = reverb(s, sr,
            echoDelay = 250, echoLevel = -20,
            reverbDelay = 70, reverbSpread = 120,
            reverbLevel = -25, reverbDensity = 50,
            filter = list(formants = NULL, lipRad = 3))
playme(s4)
spectrogram(s4, 16000, ylim = c(0, 4))

# recorded IR
s5 = soundgen(nSyl = 2, sylLen = 250, pitch = c(500, 300))
playme(s5)
spectrogram(s5, 16000, yScale = 'ERB')
s5r = reverb(s5, 16000, IR = IR)
playme(s5r)
spectrogram(s5r, 16000, yScale = 'ERB')

# add reverb to all files in a folder, save the result
reverb('~/Downloads/temp', saveAudio = TRUE)

## End(Not run)

Random draw from a truncated normal distribution

Description

rnorm_truncated generates random numbers from a normal distribution using rnorm(), but forced to remain within the specified low/high bounds. All proposals outside the boundaries (exclusive) are discarded, and the sampling is repeated until there are enough values within the specified range. Fully vectorized. Note: "truncnorm::truncnorm" is much faster, but it only accepts static low/high boundaries. rnorm_truncated2 is a simplified version of rnorm_truncated, in which values outside the bounds are simply reset to the low/high bounds. The shape of the resulting distribution is no longer Gaussian, but this is obviously much faster. Unlike in rnorm_truncated, "low" and "high" should be scalars, not vectors (ie static boundaries).

Usage

rnorm_truncated(
  n = 1,
  mean = 0,
  sd = 1,
  low = NULL,
  high = NULL,
  roundToInteger = FALSE,
  invalidArgAction = "adjust"
)

rnorm_truncated2(
  n = 1,
  mean = 0,
  sd = 1,
  low = NULL,
  high = NULL,
  roundToInteger = FALSE,
  invalidArgAction = "adjust"
)

Arguments

Value

A vector of length n.

Examples

soundgen:::rnorm_truncated(n = 3, mean = 10, sd = 5, low = 7, high = NULL,
  roundToInteger = c(TRUE, FALSE, FALSE))
soundgen:::rnorm_truncated(n = 9, mean = c(10, 50, 100), sd = c(5, 0, 20),
  roundToInteger = TRUE) # vectorized
# in case of conflicts between mean and bounds, either adjust the mean:
soundgen:::rnorm_truncated(n = 3, mean = 10, sd = .1,
  low = c(15, 0, 0), high = c(100, 100, 8), invalidArgAction = 'adjust')
#... or ignore the boundaries
soundgen:::rnorm_truncated(n = 3, mean = 10, sd = .1,
  low = c(15, 0, 0), high = c(100, 100, 8), invalidArgAction = 'ignore')

hist(soundgen:::rnorm_truncated2(n = 100, mean = 10, sd = 5, low = 7, high = NULL,
  roundToInteger = c(TRUE, FALSE, FALSE)))
hist(soundgen:::rnorm_truncated2(n = 100, mean = c(10, 50, 100), sd = c(5, 0, 20),
  roundToInteger = TRUE)) # vectorized
# in case of conflicts between mean and bounds, either sample at random
# between the boundaries...
hist(soundgen:::rnorm_truncated2(n = 100, mean = 10, sd = .1,
  low = 10, high = 15, invalidArgAction = 'adjust'))
#... or ignore the boundaries
hist(soundgen:::rnorm_truncated2(n = 100, mean = 10, sd = .1,
  low = 15, high = 100, invalidArgAction = 'ignore'))
soundgen:::rnorm_truncated2(n = 6, mean = c(0, 0, 0, 0, 0, 3),
  sd = .05, low = 0, high = 6)

Schwa-related formant conversion

Description

This function performs several conceptually related types of conversion of formant frequencies in relation to the neutral schwa sound based on the one-tube model of the vocal tract. This is useful for speaker normalization because absolute formant frequencies measured in Hz depend strongly on overall vocal tract length (VTL). For example, adult men vs. children or grizzly bears vs. dog puppies have very different formant spaces in Hz, but it is possible to define a VTL-normalized formant space that is applicable to all species and sizes. Case 1: if we know vocal tract length (VTL) but not formant frequencies, schwa() estimates formants corresponding to a neutral schwa sound in this vocal tract, assuming that it is perfectly cylindrical. Case 2: if we know the frequencies of a few lower formants, schwa() estimates the deviation of observed formant frequencies from the neutral values expected in a perfectly cylindrical vocal tract (based on the VTL as specified or as estimated from formant dispersion). Case 3: if we want to generate a sound with particular relative formant frequencies (e.g. high F1 and low F2 relative to the schwa for this vocal tract), schwa() calculates the corresponding formant frequencies in Hz. See examples below for an illustration of these three suggested uses and estimateVTL for details on the algorithm for estimating formant dispersion if VTL is not known (note that schwa calls estimateVTL with the option method = 'regression').

Usage

schwa(
  formants = NULL,
  vocalTract = NULL,
  formants_relative = NULL,
  nFormants = 8,
  interceptZero = TRUE,
  tube = c("closed-open", "open-open", "open-closed", "closed-closed"),
  speedSound = 35400,
  plot = FALSE
)

Arguments

Value

A list with the following components:

vtl_measured

VTL as provided by the user, cm

vtl_apparent

VTL estimated based on formants frequencies provided by the user, cm

dF

average distance between formants, Hz

ff_measured

formant frequencies as provided by the user, Hz

ff_schwa

formant frequencies corresponding to a neutral schwa sound in this vocal tract, Hz

ff_theoretical

formant frequencies corresponding to the user-provided relative formant frequencies, Hz

ff_relative

deviation of formant frequencies from those expected for a schwa, % (e.g. if the first ff_relative is -25, it means that F1 is 25% lower than expected for a schwa in this vocal tract)

ff_relative_semitones

deviation of formant frequencies from those expected for a schwa, semitones. Like ff_relative, this metric is invariant to vocal tract length, but the variance tends to be greater for lower vs. higher formants

ff_relative_dF

deviation of formant frequencies from those expected for a schwa, proportion of formant spacing (dF). Unlike ff_relative and ff_relative_semitones, this metric has similar variance for lower and higher formants

References

• Stevens, K. N. (2000). Acoustic phonetics (Vol. 30). MIT press.

• Anikin, A., Barreda, S. & Reby, D. (2024) A practical guide to calculating vocal tract length and scale-invariant formant patterns. Behavior Research Methods 56, 5588–5604.

See Also

estimateVTL

Examples

## CASE 1: known VTL
# If vocal tract length is known, we calculate expected formant frequencies
schwa(vocalTract = 17.5)
schwa(vocalTract = 13, nFormants = 5)
schwa(vocalTract = 13, nFormants = 5, tube = 'open-open')

## CASE 2: known (observed) formant frequencies
# Let's take formant frequencies in four vocalizations, namely
# (/a/, /i/, /mmm/, /roar/) by the same male speaker:
formants_a = c(860, 1430, 2900, NA, 5200)  # NAs are OK - here F4 is unknown
s_a = schwa(formants = formants_a, plot = TRUE)
s_a
# We get an estimate of VTL (s_a$vtl_apparent),
#   same as with estimateVTL(formants_a)
# We also get theoretical schwa formants: s_a$ff_schwa
# And we get the difference (%, semitones, or dF units) in observed vs expected
#   formant frequencies: s_a[c('ff_relative', 'ff_relative_semitones', 'ff_relative_dF')]
# [a]: F1 much higher than expected, F2 slightly lower (see plot)

formants_i = c(300, 2700, 3400, 4400, 5300, 6400)
s_i = schwa(formants = formants_i, plot = TRUE)
s_i
# The apparent VTL is slightly smaller (14.5 cm)
# [i]: very low F1, very high F2

formants_mmm = c(1200, 2000, 2800, 3800, 5400, 6400)
schwa(formants_mmm, tube = 'closed-closed', plot = TRUE)
# ~schwa, but with a closed mouth

formants_roar = c(550, 1000, 1460, 2280, 3350,
                  4300, 4900, 5800, 6900, 7900)
s_roar = schwa(formants = formants_roar, plot = TRUE)
s_roar
# Note the enormous apparent VTL (22.5 cm!)
# (lowered larynx and rounded lips exaggerate the apparent size)
# s_roar$ff_relative: high F1 and low F2-F4

schwa(formants = formants_roar[1:4], plot = TRUE)
# based on F1-F4, apparent VTL is almost 28 cm!
# Since the lowest formants are the most salient,
# the apparent size is exaggerated even further

# If you know VTL, a few lower formants are enough to get
#   a good estimate of the relative formant values:
schwa(formants = formants_roar[1:4], vocalTract = 19, plot = TRUE)
# NB: in this case theoretical and relative formants are calculated
#  based on user-provided VTL (vtl_measured) rather than vtl_apparent

## CASE 3: from relative to absolute formant frequencies
# Say we want to generate a vowel sound with F1 20% below schwa
#    and F2 40% above schwa, with VTL = 15 cm
s = schwa(formants_relative = c(-20, 40), vocalTract = 15, plot = TRUE)
# s$ff_schwa gives formant frequencies for a schwa, while
#   s$ff_theoretical gives formant frequencies for a sound with
#   target relative formant values (low F1, high F2)
schwa(formants = s$ff_theoretical)

Segment a sound

Description

Finds syllables and bursts / beats separated by background noise in long recordings (up to 1-2 hours of audio per file). Syllables are defined as continuous segments that differ from background noise based on amplitude and/or spectral contrast. Bursts are defined as local maxima in the detection contour that are high enough relative to the surrounding region. A note on long recordings: an hour of audio takes ~1 min to process, but watch your memory usage and perhaps decrease maxDur to avoid running out of RAM. This has little effect on processing speed, but dramatically cuts memory requirements. Another advantage of using shorter maxDur is that signal-noise separation may improve because noise amounts and profiles can be estimated in each chunk, adjusting for variable recording conditions over time. The downside is that small segments at chunk boundaries may be missed.

Usage

segment(
  x,
  samplingRate = NULL,
  scale = NULL,
  from = NULL,
  to = NULL,
  shortestSyl = 40,
  shortestPause = 40,
  input = c("mel", "env", "spec"),
  propNoise = NULL,
  SNR = NULL,
  specDiffMeasure = c("cosine", "logSpecDist", "specDiff"),
  amplWeight = c("add", "gate", "scale", "none"),
  noiseLevelStabWeight = c(1, 0.25),
  windowLength = 40,
  step = NULL,
  overlap = 80,
  reverb_pars = list(reverbDelay = 70, reverbSpread = 130, reverbLevel = -35,
    reverbDensity = 50, echoLevel = -Inf),
  interburst = NULL,
  peakToTrough = NULL,
  summaryFun = c("median", "sd"),
  maxDur = 120,
  ptvStep = NULL,
  ptvTime = 0.5,
  ptvFreq = 1,
  reportEvery = NULL,
  cores = 1,
  plot = FALSE,
  savePlots = FALSE,
  embed = FALSE,
  saveAudio = FALSE,
  addSilence = 50,
  main = NULL,
  xlab = "",
  ylab = NULL,
  showLegend = FALSE,
  width = 900,
  height = 500,
  units = "px",
  res = NA,
  maxPoints = c(1e+05, 5e+05),
  specPlot = list(colorTheme = "bw"),
  contourPlot = list(lty = 1, lwd = 2, col = "green"),
  sylPlot = list(lty = 1, lwd = 2, col = "blue"),
  burstPlot = list(pch = 8, cex = 3, col = "red"),
  ...
)

Arguments