2.1 Basic principles

At the heart of acoustic analysis with soundgen is the short-time Fourier transform (STFT): we look at one short segment of sound at a time (one STFT frame), analyze its spectrum using Fast Fourier Transform (FFT), and then move on to the next - perhaps overlapping - frame. As the analysis window slides along the signal, STFT shows which frequencies it contains at different points of time. The nuts and bolts of STFT are beyond the scope of this vignette, but they can be found in just about any textbook on phonetics, acoustics, digital signal processing, etc. For a quick R-friendly introduction, see seewave vignette on acoustic analysis.

Putting the spectra of all frames together, we get a spectrogram. analyze calls another function from soundgen package, spectrogram, to produce a spectrogram and then plot pitch candidates on top of it. See the examples in ?spectrogram for plot customization like color themes, contrast, brightness, etc. To analyze a sound with default settings and plot its spectrogram, all we need to specify is its sampling rate (the default in soundgen is 16000 Hz):

a1 = analyze(s1, samplingRate = 16000, plot = TRUE, ylim = c(0, 4))

## Scale not specified. Assuming that max amplitude is 1

# summary(a1)  # many acoustic predictors measured for each STFT frame
median(true_pitch)  # true value, as synthesized above

## [1] 633.2559

median(a1$pitch, na.rm = TRUE)  # our estimate

## [1] 538.42

# Pitch postprocessing is stochastic (see below), so the contour may vary.
# Many candidates are off target, mainly b/c of misleading subharmonics.

There are several key parameters that control the behavior of STFT and affect all extracted acoustic variables. The same parameters serve as arguments to spectrogram. As a result, you can immediately see what frame-by-frame input you have fed into the algorithm for acoustic analysis by visually inspecting the produced spectrogram. If you can hear f0, but can’t see individual harmonics in the spectrogram, the pitch tracker probably will not see them, either, and will therefore fail to detect f0 correctly. The first remedy is thus to adjust STFT settings, using the spectrogram for visual feedback:

• windowLength: the length of sliding STFT window. Longer windows (e.g., 40 - 50 ms) improve frequency resolution at the expense of time resolution, so they are good for detecting relatively low, slowly changing f0, as in human moans or grunts. Shorter windows (e.g., 5 - 10 ms) improve time resolution at the expense of frequency resolution, so they are good for visualizing formants or tracking high-frequency, rapidly changing f0 as in bird chirps or dolphin whistles.
• step: the step of sliding STFT window. For example, if windowLength = 50 and step = 25, each time we move the analysis frame, there is a 50% overlap with the previous frame. This introduces redundancy into the analysis, but it also - to some limited extent - improves time resolution while maintaining relatively high frequency resolution. The main cost of small steps (large overlap) is processing time, but very large overlap is not always desirable, even when processing time is not an issue. If some audio segments are problematic (e.g., very noisy), pitch contour may actually be more accurate with relatively large steps and more smoothing. It is therefore best to check the results with different steps and/or run formal optimization (remember to adjust smoothing and other postprocessing parameters together with STFT settings).
• wn: the type of windowing function used to taper the analysis frame during STFT. In practice the windowing function doesn’t seem to have a major effect on the result, as long as you choose something reasonable like gaussian, hanning, or bartlett.
• zp: zero-padding. You can use a short STFT window and improve its frequency resolution by padding each frame with zeroes. This is a computational trick that - again, to some limited extent - improves frequency resolution while maintaining relatively high time resolution.
• silence: frames with root mean square (RMS) amplitude below silence threshold are not analyzed at all. Quiet frames are harder to analyze, because their signal-to-noise ratio is lower. As a result, we want to strike a good balance. Setting silence too low (close to 0) produces a lot of garbage, as the algorithm tries to analyze frames that are essentially just background noise without any signal. Setting silence too high (close to 1) excludes too many perfectly good frames, misrepresenting the signal. In soundgen silence is dynamically updated: it can never be lower than specified, but it may be raised to the minimum root mean square amplitude of all frames, if this minimum is higher than silence. This ensures that empty frames are not analyzed in recordings with unusually high levels of steady background noise (e.g., microphone hiss).

2.2 Basic spectral descriptives

Apart from pitch tracking, analyze calculates and returns several acoustic characteristics from each non-silent STFT frame:

• time: time of the middle of each frame (ms)
• 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)
• ampl: root mean square of amplitude per frame, calculated as sqrt(mean(frame ^ 2))
• dom: lowest dominant frequency band (Hz) (see “Pitch tracking methods / Dominant frequency”)
• entropy: Weiner entropy of the spectrum of the current frame. Close to 0: pure tone or tonal sound with nearly all energy in harmonics; close to 1: white noise
• f1_freq, f1_width, …: the frequency and bandwidth of the first nFormants formants per STFT frame, as calculated by phonTools::findformants with default settings
• harmEnergy: the amount of energy in upper harmonics, namely the ratio of total spectral energy above 1.25 x f0 to the total spectral energy 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): see soundgen:::harmHeight for details
• HNR: harmonics-to-noise ratio (dB), a measure of harmonicity returned by soundgen:::getPitchAutocor() (see “Pitch tracking methods / Autocorrelation”). If HNR = 0 dB, there is as much energy in harmonics as in noise
• loudness: subjective loudness in sone, assuming a certain sound pressure level (takes into account the energy in different frequency bands as well as the sensitivity of human ears to different frequencies); see getLoudness() and section on Loudness below for details
• peakFreq: the frequency with maximum spectral energy below cutFreq (Hz)
• roughness: the amount of spectrotemporal modulation in the “roughness” zone of frequencies (estimated by modulatinSpectrum, the arguments to which are passed in roughness = list())
• quartile25, quartile50, quartile75: the 25th, 50th, and 75th quantiles of the spectrum below cutFreq (Hz) for VOICED frames
• specCentroid: the center of gravity of the frame’s spectrum below cutFreq, first spectral moment (Hz)
• specSlope: the slope of linear regression fit to the spectrum below cutFreq
• voiced: is the current STFT frame voiced? TRUE / FALSE

TIP: if voicedSeparate = TRUE, these descriptors are calculated separately for the entire sound and only for the voiced frames, resulting in extra output variables like “amplVoiced”.

2.3 Custom spectral descriptives

The function soundgen::analyze returns a few spectral descriptives that make sense for nonverbal vocalizations, but additional predictors may be useful for other applications (bird songs, non-biological sounds, etc.). One way to obtain extra predictors is to add the necessary code to the internal function soundgen:::analyzeFrame() and to soundgen::analyze(). If you want deltas, they can be extracted directly from the output of analyze(..., summaryFun = NULL). But in many cases the easiest solution may be to just extract the spectra and then process them manually, without calling analyze(). In fact, many popular spectral descriptors are mathematically trivial to derive - all you need is the spectrum for each STFT frame, or perhaps even the average spectrum of the entire sound. Here is how you can get these spectra.

For the average spectrum of an entire sound, go no further than seewave::spec or seewave::meanspec:

spec = seewave::spec(s1, f = 16000, plot = FALSE)  # FFT of the entire sound
avSpec = seewave::meanspec(s1, f = 16000, plot = FALSE)  # STFT followed by averaging
# either way, you get a dataframe with two columns: frequencies and their strength
head(avSpec)

##            x            y
## [1,] 0.00000 0.0001176617
## [2,] 0.03125 0.0001893805
## [3,] 0.06250 0.0003349440
## [4,] 0.09375 0.0009898188
## [5,] 0.12500 0.0027037798
## [6,] 0.15625 0.0038825418

If you are interested in how the spectrum changes over time, extract frame-by-frame spectra - for example, with spectrogram(..., output = 'original'):

spgm = spectrogram(s1, samplingRate = 16000, output = 'original', plot = FALSE)
# rownames give you frequencies in KHz, colnames are time stamps in ms
str(spgm)

##  num [1:400, 1:77] 2.38e-04 1.84e-04 9.22e-05 5.64e-05 3.77e-05 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:400] "0" "0.02" "0.04" "0.06" ...
##   ..$ : chr [1:77] "0" "15" "30" "45" ...

Let’s say you are working with frame-by-frame spectra and want to calculate skewness, the 66.6th percentile, and the ratio of energy above/below 500 Hz. Before you go hunting for a piece of software that returns exactly those descriptors, consider this. Once you have normalized the spectrum to add up to 1, it basically becomes a probability density function (pdf), so you can summarize it in the same way as you would any other distribution of a random variable. Look up the formulas you need and just do the raw math:

# Transform spectrum to pdf (all columns should sum to 1):
spgm_norm = apply(spgm, 2, function(x) x / sum(x))
# Set up a dataframe to store the output
out = data.frame(skew = rep(NA, ncol(spgm)),
                 quantile66 = NA,
                 ratio500 = NA)
# Process each STFT frame
for (i in 1:ncol(spgm_norm)) {
  # Absolute spectrum for this frame
  df = data.frame(
    freq = as.numeric(rownames(spgm_norm)),  # frequency (kHz)
    d = spgm_norm[, i]                       # density
  )
  # plot(df, type = 'l')
  
  # Skewness (see https://en.wikipedia.org/wiki/Central_moment)
  m = sum(df$freq * df$d)  # spectral centroid, kHz
  out$skew[i] = sum((df$freq - m)^3 * df$d)
  
  # 66.6th percentile (2/3 of density below this frequency)
  out$quantile66[i] = df$freq[min(which(cumsum(df$d) >= 2/3))]  # in kHz
  
  # Energy above/below 500 Hz
  out$ratio500[i] = sum(df$d[df$freq >= .5]) / sum(df$d[df$freq < .5])
}

## Warning in min(which(cumsum(df$d) >= 2/3)): no non-missing arguments to min;
## returning Inf

summary(out)

##       skew            quantile66        ratio500        
##  Min.   : 0.01562   Min.   :0.0200   Min.   :  0.00146  
##  1st Qu.: 0.40766   1st Qu.:0.8950   1st Qu.: 13.24145  
##  Median : 0.80829   Median :1.1200   Median : 34.35217  
##  Mean   : 1.67471   Mean   :0.9895   Mean   : 56.79107  
##  3rd Qu.: 1.20217   3rd Qu.:1.2650   3rd Qu.: 85.33104  
##  Max.   :13.68222   Max.   :1.6600   Max.   :242.96660  
##  NA's   :1          NA's   :1        NA's   :1

If you need to do this analysis repeatedly, just wrap the code into your own function that takes a wav file as input and returns all these spectral descriptives. You can also save the actual spectra of different sound files and add them up to obtain an average spectrum across multiple sound files, work with cochleograms instead of raw spectra (check out tuneR::melfcc), etc. Be your own boss!

2.4 Loudness

The digital representation of a sound is a long vector of numbers on some arbitrary scale, say [-1, 1]. Values further from zero correspond to a higher amplitude - in physical terms, to greater pertubations of sound pressure level caused by the propagating sound wave. A smoothed line following peak amplitude values is known as an amplitude envelope. However, there is no simple correspondence between the absolute height of amplitude peaks and the subjectively experienced loudness of the corresponding sound. A commonly reported measure of sound intensity is its root mean square (RMS) amplitude, which takes into account the average value of sound pressure, and not only the height of peaks. More sophisticated estimates of loudness also take into account the relative sensitivity of human hearing to different frequencies, masking of adjacent tones in the time and frequency domains, etc.

To illustrate the differences between these estimates, let’s look at a pure tone sweeping with fixed absolute amplitude from 100 to 4000 Hz over 2 s:

dur = 2  # 2 s duration
samplingRate = 16000
f0 = seq(100, 8000, length.out = samplingRate * dur)
sweep = sin(2 * pi * cumsum(f0) / samplingRate)
# playme(sweep)
# spectrogram(sweep, 16000)
# plot(sweep, type = 'l')

Smoothed absolute amplitude envelope (flat):

seewave::env(sweep, f = samplingRate, envt = 'abs', msmooth=c(50, 0))

RMS amplitude per STFT frame, as returned by analyze(), column “ampl”:

a = analyze(sweep, samplingRate = samplingRate, pitchMethods = NULL, plot = FALSE)

## Scale not specified. Assuming that max amplitude is 1

plot(seq(0, dur, length.out = length(a$ampl)), a$ampl, type = 'b', xlab = 'Time, s')

An estimate of subjectively experienced loudness in sone, column “loudness”:

plot(seq(0, dur, length.out = length(a$loudness)), a$loudness, type = 'b', xlab= 'Time, s')

Soundgen also has a dedicated function for calculating the loudness and plotting the output, getLoudness(). Loudness values are overlaid on the spectrogram - observe how the loudness peaks as f0 reaches about 2-3 kHz and then drops. The absolute values in sone are only an approximation, since they are dictated by the playback device (e.g. your headphones), but the change of loudness within one sound, or across different sounds analyzed with the same settings, is informative.

l = getLoudness(sweep, samplingRate = samplingRate)

## Warning in getLoudness(sweep, samplingRate = samplingRate): Scale not specified.
## Assuming that max amplitude is 1

2.5 Pitch tracking

If you look at the source code of soundgen::analyze() and embedded functions, you will see that almost all of this code deals with a single acoustic characteristic: fundamental frequency (f0) or its perceptual equivalent, pitch. That’s because pitch is both highly salient to listeners and notoriously difficult to measure accurately. The approach followed by soundgen’s pitch tracker is to use several different estimates of f0, each of which is better suited to certain types of sounds. You can use any pitch tracker individually, but their output is also automatically integrated and postprocessed so as to generate the best overall estimate of frame-by-frame pitch. There are four currently implemented classes of pitch estimates in soundgen: autocorrelation, lowest dominant frequency, cepstrum, and spectrum (ratios of harmonics). These four methods of pitch estimation are not treated as completely independent in soundgen. Autocorrelation is performed first to provide an initial guess at the likely pitch and harmonics-to-noise ratio (HNR) of an STFT frame, and then this information is used to adjust the expectations of the cepstral and spectral algorithms. In particular, if autocorrelation suggests that the pitch is high, confidence in cepstral estimates is attenuated; and if autocorrelation suggests that HNR is low, thresholds for spectral peak detection are raised, making spectral pitch estimates more conservative.

The plot below shows a spectrogram of the sound with overlaid pitch candidates generated by five different methods (listed in pitchMethods), with a very vague prior - that is, with no specific expectations regarding the true range of pitch values. The size of each point shows the certainty of estimation: smaller points are calculated with lower certainty and have less weight when all candidates are integrated into the final pitch contour (blue line).

a = analyze(s1, samplingRate = 16000, priorSD = 24, ylim = c(0, 4),
            pitchMethods = c('autocor', 'cep', 'dom', 'spec', 'hps'))

## Scale not specified. Assuming that max amplitude is 1

Different pitch tracking methods have their own pros and cons. Cepstrum is helpful for speech but pretty useless for high-frequency whistles or screams, harmonic product spectrum (hps) is easily mislead by subharmonics (as in this example), lowest dominant frequency band (dom) can’t handle low-frequency wind noise, etc. The default is to use “dom” and “autocor” as the most generally applicable, but you can experiment with all methods and check which ones perform best with the specific type of audio that you are analyzing. Each method can also be fine-tuned (see below), but first it is worth considering the general pitch-related settings.

par(mfrow = c(1, 2))
# default prior in soundgen
getPrior(priorMean = 300, priorSD = 6)
# narrow peak at 2 kHz
getPrior(priorMean = 2000, priorSD = 1)
par(mfrow = c(1, 1))

a = analyze(s1, samplingRate = 16000, 
            plot = TRUE, ylim = c(0, 4), priorMean = NA,
            pitchMethods = 'autocor',
            pitchAutocor = list(autocorThres = .45, 
                                # + plot pars if needed
                                col = 'green'),
            nCands = 3)

## Scale not specified. Assuming that max amplitude is 1

a = analyze(s1, 
            samplingRate = 16000, ylim = c(0, 4), priorMean = NA,
            pitchMethods = 'dom',
            pitchDom = list(domThres = .1, domSmooth = 500, cex = 1.5))

## Scale not specified. Assuming that max amplitude is 1

a = analyze(s1, 
            samplingRate = 16000, ylim = c(0, 4), priorMean = NA,
            pitchMethods = 'cep',
            pitchCep = list(cepThres = .3),
            nCands = 2)

## Scale not specified. Assuming that max amplitude is 1

a = analyze(s1, 
            samplingRate = 16000, plot = TRUE, ylim = c(0, 4), priorMean = NA,
            pitchMethods = 'spec',
            pitchSpec = list(specThres = .2, specPeak = .1, cex = 2),
            nCands = 2)

## Scale not specified. Assuming that max amplitude is 1

a = analyze(s1, 
            samplingRate = 16000, plot = TRUE, ylim = c(0, 4), priorMean = NA,
            pitchMethods = 'hps',
            pitchHps = list(hpsNum = 2, # try 8 or so to measure subharmonics
                            hpsThres = .2))

## Scale not specified. Assuming that max amplitude is 1

s_withf0 = soundgen(sylLen = 600, pitch = 300,
              rolloffExact = c(1, 1, 1, 1), formants = NULL, lipRad = 0)
# playme(s_withf0)
seewave::meanspec(s_withf0, f = 16000, dB = 'max0', flim = c(0, 3))

s_withoutf0 = soundgen(sylLen = 600, pitch = 300,
              rolloffExact = c(0, 1, 1, 1), formants = NULL, lipRad = 0)
# playme(s_withoutf0)  # you can clearly hear the difference
seewave::meanspec(s_withoutf0, f = 16000, dB = 'max0', flim = c(0, 3))

a_withoutf0 = analyze(s_withoutf0, 16000, 
             pitchMethods = c('autocor', 'dom', 'cep', 'spec', 'hps'),
             ylim = c(0, 2), dynamicRange = 60, osc = FALSE, priorMean = NA)

## Scale not specified. Assuming that max amplitude is 1

2.6 Postprocessing of pitch contour

Pitch postprocessing in soundgen includes a whole battery of distinct operations through which the pitch candidates generated by one or more tracking methods are integrated into the final pitch contour. We will look at them one by one, in the order in which they are performed in analyze. But first of all, here is how to disable them all:

a = analyze(
  s1, 
  samplingRate = 16000, plot = TRUE, ylim = c(0, 4), priorMean = NA,
  shortestSyl = 0,        # any length of voiced fragments
  interpolWin = 0,        # don't interpolate missing f0 values
  pathfinding = 'none',   # don't look for optimal path through candidates
  snakeStep = 0,          # don't run the snake
  smooth = 0              # don't run median smoothing
)       

## Scale not specified. Assuming that max amplitude is 1

When the sound is not too tricky and enough pitch candidates are available, postprocessing actually makes little difference. In terms of the accuracy of median estimate of f0, you are likely to get a good result even with postprocessing is completely disabled. However, if you are interested in the actual intonation contours, not just the global average, postprocessing can help a lot.

a1 = analyze(s1, samplingRate = 16000, priorMean = NA,
             pitchMethods = 'cep', pitchCep = list(cepThres = .4), step = 25,
             snakeStep = 0, smooth = 0,
             interpolWin = 0,   # disable interpolation
             pathfinding = 'none',  
             summaryFun = NULL,
             plot = FALSE)

## Scale not specified. Assuming that max amplitude is 1

a2 = analyze(s1, samplingRate = 16000, priorMean = NA,
             pitchMethods = 'cep', pitchCep = list(cepThres = .4), step = 25,
             pathfinding = 'none',
             snakeStep = 0, smooth = 0, 
             summaryFun = NULL,
             plot = FALSE)  

## Scale not specified. Assuming that max amplitude is 1

plot(a1$time, a1$pitch, type = 'l', xlab = 'Time, ms', ylab = 'Pitch, Hz')
points(a2$time, a2$pitch, type = 'l', col = 'red', lty = 3)

a1 = analyze(s1, samplingRate = 16000, priorMean = NA,
             pitchMethods = 'cep', pitchCep = list(cepThres = .15), nCands = 3,
             snakeStep = 0, smooth = 0, interpolTol = Inf,
             certWeight = 0,  # minimize pitch jumps
             main = 'Minimize jumps', 
             showLegend = FALSE, osc = FALSE, ylim = c(0, 3))  

## Scale not specified. Assuming that max amplitude is 1

a2 = analyze(s1, samplingRate = 16000, priorMean = NA,
             pitchMethods = 'cep', pitchCep = list(cepThres = .15), nCands = 3,
             snakeStep = 0, smooth = 0, interpolTol = Inf,
             certWeight = 1,  # minimize deviation from high-certainty candidates
             main = 'Pass through top cand-s', 
             showLegend = FALSE, osc = FALSE, ylim = c(0, 3))

## Scale not specified. Assuming that max amplitude is 1

a1 = analyze(s1, samplingRate = 16000, plot = FALSE, priorMean = NA,
             pitchMethods = 'cep', pitchCep = list(cepThres = .2), nCands = 2,
             pathfinding = 'none', smooth = 0, interpolTol = Inf,
             certWeight = 0.1,  # like pathfinding, the snake is affected by certWeight
             snakeStep = 0.05, snakePlot = TRUE)

## Scale not specified. Assuming that max amplitude is 1