Measuring schroeder periods

Fine scale acoustic perception in zebra finches

Source code and data found at https://github.com/maRce10/acoustic-fine-features-zebra-finch

 

1 Purpose

• Detect start and end of periods in a schroeder using amplitude cross-correlation
• Compare methods for measuring fine scale structural variation and periodicity in schroeders

 

2 Report overview

• Synthetizing schroeders
• Detecting periodicity
• Measuring schroeder dissimilarity
• Takeaways

 

source("~/Dropbox/Projects/geographic_call_variation_yellow-naped_amazon/scripts/MRM2.R")

warbleR_options(parallel = 22)

wave_col <- viridis(10)[7]

3 Synthetizing schroeders

3.1 Function to make schroeders

qwindpc <- function(rftime, srate, sig) {
    ramppts <- round((rftime/1000) * srate)
    hold1 <- seq(0, 1, length.out = ramppts + 1)
    onramp <- sin(0.5 * pi * hold1)^2
    offramp <- cos(0.5 * pi * hold1)^2
    steady <- rep(1, length(sig) - (2 * ramppts) - 2)
    wind <- c(onramp, steady, offramp)

    if (length(wind) > length(sig)) {
        print("Window not applied, block too short for window")
    } else {
        sig <- wind * sig
    }

    return(sig)
}


make_schroeder <- function(samp.rate = 44100, f0, dur = 1000, save.wave = FALSE,
    plot = TRUE, n.components, scalar = 1, color = wave_col, path = ".",
    file.name = NULL, random.start = FALSE, mfrow = c(2, 1), par = TRUE,
    cex.oscillo = 1, ...) {

    # duplicate duration if starting randomly
    dur.mult <- if (random.start)
        2 else 1

    # set time instances at which samples will be 'taken'
    t <- seq(1/samp.rate, 1, by = 1/samp.rate)
    sumwavepos <- rep(0, samp.rate)

    durpts <- round((dur * dur.mult/1000) * samp.rate)

    phase <- rep(0, n.components)
    amplin <- rep(1, n.components)
    f <- rep(0, n.components)

    Ns <- 1:n.components

    for (n in Ns) {
        compnum <- n
        phase[compnum] <- scalar * pi * compnum * (compnum - 1)/n.components
        f[n] <- f0 + ((compnum - 1) * f0)
    }

    wave <- matrix(0, nrow = n.components, ncol = length(t))

    for (i in Ns) {
        wave[i, ] <- amplin[i] * cos((2 * pi * f[i] * t) + phase[i])
        sumwavepos <- sumwavepos + wave[i, ]
    }

    # make sumwapos longer so it fit the target duration
    org.sumwavepos <- sumwavepos

    while (length(sumwavepos) < durpts) sumwavepos <- c(sumwavepos,
        org.sumwavepos)

    posschr <- sumwavepos[1:durpts]/max(abs(sumwavepos[1:durpts]))

    # plot(t[1:length(posschr)], posschr, type = 'l')

    spectr <- fft(sumwavepos)
    power <- abs(spectr)^2
    sample <- length(t)/max(t)
    freq <- (1:length(t)) * sample/length(t)
    powerdb <- 10 * log10(power[1:(length(t)/2)])

    if (random.start) {
        strt <- sample(1:(length(posschr)/2), 1)
        posschr <- posschr[strt:(strt + (length(posschr)/2) - 1)]
    }

    posschr <- qwindpc(10, samp.rate, posschr)
    posschr <- posschr/max(abs(posschr))

    wave_obj <- tuneR::normalize(tuneR::Wave(posschr, samp.rate = samp.rate,
        bit = 16), unit = "16")

    if (plot) {
        opar <- par()
        opar$cin <- opar$cra <- opar$sci <- opar$cxy <- opar$din <- opar$csi <- opar$page <- NULL

        # on.exit(par(opar))
        if (par)
            par(mfrow = mfrow, mar = c(5, 4, 0, 0) + 0.1)

        plot(freq[1:(length(t)/2)], powerdb, type = "l", xlim = c(0,
            5200), xlab = "Frequency (Hz)", ylab = "Power(dB)", col = color,
            ...)

        seewave::oscillo(wave = wave_obj, colwave = color, cexlab = cex.oscillo,
            cexaxis = cex.oscillo)
    }

    # Save audio to a file
    if (save.wave) {
        if (is.null(file.name)) {
            file.name <- paste0("f0-", f0, "_ncomp-", n.components,
                "_", if (scalar == 1)
                  "pos" else "neg", ".wav")
        }
        writeWave(object = wave_obj, filename = file.path(path, file.name),
            extensible = FALSE)

    } else {
        return(wave_obj)
    }
}

3.1.1 Schroeder at 700 Hz with 7 components (harmonics)

Positive

ms <- make_schroeder(n.components = 7, samp.rate = 44100, f0 = 400,
    dur = 25, save.wave = FALSE, plot = TRUE, scalar = 1, color = viridis(10)[3],
    random.start = FALSE)

Negative

ms <- make_schroeder(n.components = 7, samp.rate = 44100, f0 = 400,
    dur = 25, save.wave = FALSE, plot = TRUE, scalar = -1, color = viridis(10)[3])

3.1.2 Schroeder with 10 components

Positive

ms <- make_schroeder(n.components = 10, samp.rate = 44100, f0 = 400,
    dur = 25, save.wave = FALSE, plot = TRUE, scalar = 1, color = viridis(10)[3])

Negative

ms <- make_schroeder(n.components = 10, samp.rate = 44100, f0 = 400,
    dur = 25, save.wave = FALSE, plot = TRUE, scalar = -1, color = viridis(10)[3])

Randomly starting at any point

ms <- make_schroeder(n.components = 10, samp.rate = 44100, f0 = 400,
    dur = 25, save.wave = FALSE, plot = TRUE, scalar = -1, color = viridis(10)[3],
    random.start = TRUE)

Make example video (or gif):

ani.options(ani.res = 200, ani.width = 1700, ani.height = 1000, ani.type = "png")

## make sure ImageMagick has been installed in your system gifs
## <- function(lng_dat, type = 'gif'){

# anni_fun <- if(type == 'gif') saveGIF else saveVideo


saveGIF(expr = {

    for (i in 1:13) {
        par(mfcol = c(2, 2), mar = c(5, 5, 3, 2) + 0.1)

        ms <- make_schroeder(n.components = i, samp.rate = 44100,
            f0 = 400, dur = 25, save.wave = FALSE, plot = TRUE, scalar = 1,
            color = viridis(10)[3], random.start = FALSE, par = FALSE,
            main = paste0("Positive, ", i, " components"), cex.main = 2,
            cex.lab = 2, cex.oscillo = 2)

        ms <- make_schroeder(n.components = i, samp.rate = 44100,
            f0 = 400, dur = 25, save.wave = FALSE, plot = TRUE, scalar = -1,
            color = viridis(10)[8], random.start = FALSE, par = FALSE,
            main = paste0("Negative, ", i, " components"), cex.main = 2,
            cex.lab = 2, cex.oscillo = 2)

    }


}, video.name = file.path("./output/movies", "schroeders_by_sign_and_harmonics2"),
    interval = 1)
# }, movie.name = paste0(z, '.', type), interval = 0.05) }

schroeders_by_sign_and_harmonics

3.2 Make schroeders

Full factorial design of the following features:

• Components: from 5 to 25 (21 values)
• Fundamental frequency: from 200 to 800 (7 values)
• Phase: 1 and -1
• 294 schroeders

Create schroeders as individual sound files:

• Fixed start

shr_tab <- expand.grid(type = c(-1, 1), components = seq(5, 25, 1),
    f0 = seq(200, 800, 100))

nrow(shr_tab)


shr_tab$name <- paste0(f_pad_zero(x = 1:nrow(shr_tab), width = 3),
    "_f0-", shr_tab$f0, "_ncomp-", shr_tab$components, "_", ifelse(shr_tab$type ==
        1, "pos", "neg"), ".wav")

out <- warbleR:::pblapply_wrblr_int(seq_len(nrow(shr_tab)), function(x) {
    make_schroeder(n.components = shr_tab$components[x], samp.rate = 44100,
        f0 = shr_tab$f0[x], dur = 1000, save.wave = TRUE, plot = FALSE,
        scalar = shr_tab$type[x], color = "red", path = "./data/processed/schroeders",
        file.name = shr_tab$name[x])
})

• Random start

shr_tab <- expand.grid(type = c(-1, 1), components = seq(5, 25, 1),
    f0 = seq(200, 800, 100))

nrow(shr_tab)


shr_tab$name <- paste0(f_pad_zero(x = 1:nrow(shr_tab), width = 3),
    "_f0-", shr_tab$f0, "_ncomp-", shr_tab$components, "_", ifelse(shr_tab$type ==
        1, "pos", "neg"), ".wav")

out <- warbleR:::pblapply_wrblr_int(seq_len(nrow(shr_tab)), function(x) {
    make_schroeder(n.components = shr_tab$components[x], samp.rate = 44100,
        f0 = shr_tab$f0[x], dur = 1000, save.wave = TRUE, plot = FALSE,
        scalar = shr_tab$type[x], color = "red", path = "./data/processed/random_start_schroeders",
        file.name = shr_tab$name[x], random.start = TRUE)
})

Put Schroeders into a single sound file

# fix start
est_schr <- warbleR::selection_table(whole.recs = TRUE, path = "./data/processed/schroeders",
    extended = TRUE, confirm.extended = FALSE)

est_schr$bottom.freq <- (as.numeric(gsub("f0-", "", (sapply(strsplit(est_schr$sound.files,
    split = "_"), "[[", 2)))) - 50)/1000

est_schr$top.freq <- freq_range(X = est_schr, threshold = 2, fsmooth = 0.8,
    parallel = 13)$top.freq

est_schr$f0 <- gsub("f0-", "", sapply(strsplit(est_schr$sound.files,
    split = "_"), "[[", 2))

est_schr$label <- gsub("os.wav_1|eg.wav_1", "", sapply(strsplit(est_schr$sound.files,
    split = "ncomp-"), "[[", 2))

saveRDS(est_schr, "./data/processed/extended_selection_table_schroeders.RDS")

master <- master_sound_file(X = est_schr, file.name = "schroeder_master",
    dest.path = "./data/processed/", gap.duration = 0.8, cex = 14)

write.csv(master, "./data/processed/master_annotations_schroeders.csv",
    row.names = FALSE)

master$f0 <- c("start_marker", paste("f0", est_schr$f0), "")


master$label <- c("marker", est_schr$label, "marker")

Rraven::exp_raven(master, path = "./data/processed/", sound.file.path = "./data/processed/")

warbleR::full_spectrograms(master, path = "./data/processed/", sxrow = duration(readWave("./data/processed/schroeder_master.wav"))/12,
    rows = 12, horizontal = TRUE, dest.path = "./data/processed/spectrograms",
    collevels = seq(-150, 0, 5), labels = "label", song = "f0", wl = 2000)


# random start
est_schr <- warbleR::selection_table(whole.recs = TRUE, path = "./data/processed/random_start_schroeders",
    extended = TRUE, confirm.extended = FALSE)

est_schr$bottom.freq <- (as.numeric(gsub("f0-", "", (sapply(strsplit(est_schr$sound.files,
    split = "_"), "[[", 2)))) - 50)/1000

est_schr$top.freq <- freq_range(X = est_schr, threshold = 2, fsmooth = 0.8,
    parallel = 13)$top.freq

est_schr$f0 <- gsub("f0-", "", sapply(strsplit(est_schr$sound.files,
    split = "_"), "[[", 2))

est_schr$label <- gsub("os.wav_1|eg.wav_1", "", sapply(strsplit(est_schr$sound.files,
    split = "ncomp-"), "[[", 2))

saveRDS(est_schr, "./data/processed/extended_selection_table_random_start_schroeders.RDS")

master <- master_sound_file(X = est_schr, file.name = "random_start_schroeder_master",
    dest.path = "./data/processed/", gap.duration = 0.8, cex = 14)

write.csv(master, "./data/processed/master_annotations_random_start_schroeders.csv",
    row.names = FALSE)

master$f0 <- c("start_marker", paste("f0", est_schr$f0), "")


master$label <- c("marker", est_schr$label, "marker")

Rraven::exp_raven(master, file.name = "random_start_schroeders", path = "./data/processed/",
    sound.file.path = "./data/processed/")

master sound file

4 Detecting periodicity

• Create a function that returns each detected segment in a list
• Two methods:
  • Empirical Mode Decomposition (EMD)
  • Time Autocorrelation (ac)
• The function can plot the mean period +/- standard deviation

est_schr <- readRDS("./data/processed/extended_selection_table_schroeders.RDS")

source("~/Dropbox/Projects/acoustic_fine_features_zebra_finch/scripts/mean_segment.R")

exp_est(est_schr, path = "./data/raw/extended_selection_table_schroeders_clips",
    file.name = "annotations")

est_schr <- resample_est(est_schr, bit.depth = 32)

# add silence
for (i in seq_len(nrow(est_schr))) {
    wv <- read_sound_file(est_schr, index = i)
    wv <- normalize(wv, unit = "32")

    # make it a 200 ms cut
    wv <- cutw(wave = wv, from = 0, to = 0.1, output = "Wave")
    # sil <- runif(min = min(wv@left) / 7, max = max(wv@left) /
    # 7, n = wv@samp.rate / 10) wv@left <- c(sil, wv@left, sil)
    # wv <- addsilw(wave = wv, at = 'end', d = 0.1, f = 44100,
    # output = 'Wave', ) wv <- addsilw(wave = wv, at = 'start',
    # d = 0.1, f = 44100, output = 'Wave')
    spectro(wv)

    filename <- file.path("./data/raw/extended_selection_table_schroeders_clips/100ms",
        gsub("_1$", "", est_schr$sound.files[i]))

    writeWave(object = wv, filename = filename, extensible = FALSE)
}

# fix_wavs(path =
# './data/raw/extended_selection_table_schroeders_clips/200ms')

4.1 Test function

mean period +/- standard deviation using autocorrelation

est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")

wave <- read_wave(est_schr, index = 40)

ms <- mean_segment(wave, thinning = 1, col = wave_col, mean = FALSE,
    plot = TRUE)

4.2 Get mean schroeder cycles

4.2.1 EMD method

est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")

mean_schroeders <- warbleR:::pblapply_wrblr_int(cl = 20, seq_len(nrow(est_schr)),
    function(x) {
        wave <- read_wave(est_schr, index = x)
        seg <- try_na(mean_segment(wave, plot = FALSE, mean = FALSE,
            type = "EMD"))

        return(seg)
    })

names(mean_schroeders) <- paste(est_schr$f0, est_schr$label, sep = "-")

saveRDS(mean_schroeders, "./data/processed/mean_random_start_schroeders_emd.RDS")

mean_schroeders <- readRDS("./data/processed/mean_random_start_schroeders_emd.RDS")

mean_schroeders <- mean_schroeders[!sapply(mean_schroeders, function(x) is.na(x[[1]][1]))]

mean_schroeders_list <- lapply(seq_len(length(mean_schroeders)), function(x) {
    data.frame(schroeder = names(mean_schroeders)[x], time = seq(0,
        1, length.out = nrow(mean_schroeders[[x]])), mean.amp = rowMeans(mean_schroeders[[x]]),
        sd.amp = apply(mean_schroeders[[x]], 1, sd))
})

mean_schroeders_df <- do.call(rbind, mean_schroeders_list)

ggplot(data = mean_schroeders_df[mean_schroeders_df$schroeder %in%
    unique(mean_schroeders_df$schroeder)[1:40], ], mapping = aes(x = time,
    y = mean.amp)) + geom_line(color = wave_col) + geom_ribbon(aes(ymin = mean.amp -
    sd.amp, ymax = mean.amp + sd.amp), alpha = 0.2) + theme_classic(base_size = 5) +
    facet_wrap("~ schroeder", ncol = 5, scales = "free_y")

4.2.2 Autocorrelation method

est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")

mean_schroeders <- warbleR:::pblapply_wrblr_int(cl = 20, seq_len(nrow(est_schr)),
    function(x) {
        wave <- read_wave(est_schr, index = x)
        seg <- try_na(mean_segment(wave, plot = FALSE, mean = FALSE,
            type = "ac", thinning = 0.8))

        return(seg)
    })

names(mean_schroeders) <- paste(est_schr$f0, est_schr$label, sep = "-")

saveRDS(mean_schroeders, "./data/processed/mean_random_start_schroeders_ac.RDS")

mean_schroeders <- readRDS("./data/processed/mean_random_start_schroeders_ac.RDS")

mean_schroeders <- mean_schroeders[!sapply(mean_schroeders, function(x) is.na(x[[1]][1]))]

mean_schroeders_list <- lapply(seq_len(length(mean_schroeders)), function(x) {
    data.frame(schroeder = names(mean_schroeders)[x], time = seq(0,
        1, length.out = nrow(mean_schroeders[[x]])), mean.amp = rowMeans(mean_schroeders[[x]]),
        sd.amp = apply(mean_schroeders[[x]], 1, sd))
})

mean_schroeders_df <- do.call(rbind, mean_schroeders_list)

ggplot(data = mean_schroeders_df[mean_schroeders_df$schroeder %in%
    unique(mean_schroeders_df$schroeder)[1:40], ], mapping = aes(x = time,
    y = mean.amp)) + geom_line(color = wave_col) + geom_ribbon(aes(ymin = mean.amp -
    sd.amp, ymax = mean.amp + sd.amp), alpha = 0.2) + theme_classic(base_size = 5) +
    facet_wrap("~ schroeder", ncol = 5, scales = "free_y")

5 Measuring schroeder dissimilarity

5.1 Dynamic-time warping pairwise distance

• Both Schroeders have the same length
• One is duplicated and the other one is slide across the duplicated one
• The minimum DTW distance is kept as a dissimilarity measure

mean_schroeders <- readRDS("./data/processed/mean_random_start_schroeders_ac.RDS")

nms <- names(mean_schroeders)
# nms <- grep(pattern = '200|400|600|800', nms, value = TRUE)

cmbs <- t(combn(nms, 2))


min_dist_l <- pbapply::pbsapply(cl = 22, 1:nrow(cmbs), function(x) {
    s1 <- rowMeans(mean_schroeders[[cmbs[x, 1]]])
    s2 <- rowMeans(mean_schroeders[[cmbs[x, 2]]])

    # make same length if (length(s1) != length(s2))
    s1 <- approx(s1, n = 100)$y
    s2 <- approx(s2, n = 100)$y

    # duplicate 1
    s1 <- rep(s1, 2)

    # run dtw over longer vector
    dists <- vapply(seq_len(length(s1) - length(s2)), function(x) {
        segment <- s1[x:min(c(x + length(s2) - 1), length(s1))]

        dtw_dist <- warbleR::try_na(dtw::dtwDist(mx = rbind(s2, segment)))
        return(dtw_dist[1, 2])
    }, FUN.VALUE = numeric(1))

    return(data.frame(schr1 = cmbs[x, 1], schdr2 = cmbs[x, 2], min(dists)))
})

min_dists <- do.call(rbind, min_dist_l)

min_dists <- as.data.frame(matrix(min_dists[, 1], ncol = 3, byrow = TRUE))

names(min_dists) <- c("schr1", "schr2", "dist")

min_dists$dist <- as.numeric(min_dists$dist)

saveRDS(min_dists, "./data/processed/dtw_distance_ac_segments_random_start.RDS")

6 Compare dissimilarity between schroeders using different methods

We try 11 methods for measuring acoustic structure:

• Dynamic time warping (DTW) distances from schroeder waveform
• Raven’s spectrographic features
• warbleR’s spectrographic features
• warbleR’s MFCC features
• warbleR’s cross-correlation of fourier spectrograms
• Raven’s cross-correlation of fourier spectrograms
• Raven’s waveform correlation
• warbleR’s cross-correlation of mel-frequency spectrograms
• DTW on warbleR’s dominant frequency contours
• DTW on Raven’s peak frequency contours
• Sound Analysis Pro similarity

Statistical modeling:

• Multiple Regression on distance Matrices
• Model:
  \[\begin{align*} Dissimilarity &\sim frequency + components + phase \end{align*}\]
• Response values scaled to make effect sizes comparable across models
• Predictors were coded as pairwise binary matrices in which 0 means that calls in a dyad belong to the same level and 1 means calls belong to different levels

6.1 Waveform DTW

min_dists <- readRDS("./data/processed/dtw_distance_ac_segments_random_start.RDS")

min_dists$dist <- scale(min_dists$dist)

dist_tri <- PhenotypeSpace::rectangular_to_triangular(min_dists)

freq_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
    "-"), "[[", 1)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = gsub("_n|_p",
    "", sapply(strsplit(rownames(dist_tri), "-"), "[[", 2))))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
    "_"), "[[", 2)))

rect_var <- cbind(as.dist(dist_tri), freq_bi_tri, comp_bi_tri, sign_bi_tri)

colnames(rect_var) <- c("dtw_dist", "frequency", "components", "sign")

dtw_wv_mod <- MRM2(formula = dtw_dist ~ frequency + components + sign,
    nperm = 10000, mat = rect_var)

saveRDS(dtw_wv_mod, "./data/processed/matrix_correlation_dtw_distance_random_start.RDS")

(dtw_wv_mod <- readRDS("./data/processed/matrix_correlation_dtw_distance_random_start.RDS"))

$coef
             dtw_dist  pval
Int        -1.5156699 1e-04
frequency   0.1859515 1e-04
components  1.0638637 1e-04
sign        0.6755587 1e-04

$r.squared
       R2      pval 
0.1624519 0.0001000 

$F.test
        F    F.pval 
2803.4517    0.0001 

$n
[1] 295

6.2 Raven’s spectrographic features

rav_dat <- imp_raven(path = "./data/processed", files = "random_start_schroeder_master_raven_measurements.txt",
    warbler.format = TRUE, all.data = TRUE)

rav_dat <- rav_dat[grep(pattern = "marker", rav_dat$f0, invert = TRUE),
    ]

rav_dat <- rav_dat[, grep("contour|PFC Slope", ignore.case = TRUE,
    x = names(rav_dat), invert = TRUE)]

pca <- prcomp(x = rav_dat[, names(rav_dat) %in% c("Agg Entropy (bits)",
    "Avg Entropy (bits)", "Avg Power Density (dB FS/Hz)", "BW 50% (Hz)",
    "BW 90% (Hz)", "Center Freq (Hz)", "Center Time Rel.", "Delta Freq (Hz)",
    "Dur 50% (s)", "Dur 90% (s)", "Energy (dB FS)", "Freq 25% (Hz)",
    "Freq 5% (Hz)", "Freq 75% (Hz)", "Freq 95% (Hz)", "Inband Power (dB FS)",
    "Max Entropy (bits)", "Max Freq (Hz)", "Min Entropy (bits)", "Peak Freq (Hz)",
    "PFC Avg Slope (Hz/ms)", "PFC Max Freq (Hz)", "PFC Max Slope (Hz/ms)",
    "PFC Min Freq (Hz)", "PFC Min Slope (Hz/ms)", "PFC Num Inf Pts",
    "Peak Power Density (dB FS/Hz)", "Peak Time (s)", "Peak Time Relative",
    "Time 25% (s)", "Time 25% Rel.", "Time 5% (s)", "Time 5% Rel.",
    "Time 75% (s)", "Time 75% Rel.", "Time 95% (s)", "Time 95% Rel.")],
    scale. = TRUE)

rav_dat$pc1 <- pca$x[, 1]
rav_dat$pc2 <- pca$x[, 2]
rav_dat$comp <- sapply(strsplit(rav_dat$label, "_"), "[[", 1)
rav_dat$sign <- sapply(strsplit(rav_dat$label, "_"), "[[", 2)
rav_dat$label <- paste(rav_dat$f0, rav_dat$label, sep = "-")

dist_tri <- dist(rav_dat[, c("pc1", "pc2")])

freq_bi_tri <- as.dist(binary_triangular_matrix(group = rav_dat$f0))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = rav_dat$comp))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = rav_dat$sign))

rect_var <- cbind(as.dist(scale(as.matrix(dist_tri))), freq_bi_tri,
    comp_bi_tri, sign_bi_tri)

colnames(rect_var) <- c("rav_dist", "frequency", "components", "sign")


rav_mod <- MRM2(formula = rav_dist ~ frequency + components + sign,
    nperm = 10000, mat = rect_var)

saveRDS(rav_mod, "./data/processed/matrix_correlation_raven_measurements_distance_random_start.RDS")

(rav_mod <- readRDS("./data/processed/matrix_correlation_raven_measurements_distance_random_start.RDS"))

$coef
               rav_dist   pval
Int        -1.308738783 0.0001
frequency   0.879062058 0.0001
components  0.655162497 0.0001
sign       -0.004605695 0.5591

$r.squared
        R2       pval 
0.09958528 0.00010000 

$F.test
        F    F.pval 
1587.7273    0.0001 

6.3 warbleR’s spectrographic features

est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")

wsp <- spectro_analysis(est_schr, harmonicity = TRUE, nharmonics = 5,
    parallel = 20)

# keep columns with no NAs
wsp <- wsp[complete.cases(wsp), ]

pca <- prcomp(x = wsp[, -c(1:3)])
wsp$pc1 <- pca$x[, 1]
wsp$pc2 <- pca$x[, 2]

wsp_tri <- dist(wsp[, c("pc1", "pc2")])

freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(wsp$sound.files,
    "_"), "[[", 2), start = 4, 6)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(wsp$sound.files,
    "_"), "[[", 3), start = 7, 8)))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
    wsp$sound.files), "pos", "neg")))

rect_var <- cbind(as.dist(scale(as.matrix(wsp_tri))), freq_bi_tri,
    comp_bi_tri, sign_bi_tri)

colnames(rect_var) <- c("wrbl_sp", "frequency", "components", "sign")

wrbl_mod <- MRM2(formula = wrbl_sp ~ frequency + components + sign,
    nperm = 10000, mat = rect_var)

saveRDS(wrbl_mod, "./data/processed/matrix_correlation_warbler_measurements_distance_random_start.RDS")

(wrbl_mod <- readRDS("./data/processed/matrix_correlation_warbler_measurements_distance_random_start.RDS"))

$coef
                wrbl_sp   pval
Int        -1.310541864 0.0001
frequency   1.536291158 0.0001
components  0.104840203 0.0002
sign       -0.005210514 0.5868

$r.squared
       R2      pval 
0.3557964 0.0001000 

$F.test
        F    F.pval 
4116.8767    0.0001 

6.4 warbleR’s MFCC features

est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")

wsp <- mfcc_stats(est_schr)

# keep columns with no NAs
wsp <- wsp[complete.cases(wsp), ]

pca <- prcomp(x = wsp[, -c(1:3)])
wsp$pc1 <- pca$x[, 1]
wsp$pc2 <- pca$x[, 2]

wsp_tri <- dist(wsp[, c("pc1", "pc2")])

freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(wsp$sound.files,
    "_"), "[[", 2), start = 4, 6)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(wsp$sound.files,
    "_"), "[[", 3), start = 7, 8)))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
    wsp$sound.files), "pos", "neg")))

rect_var <- cbind(as.dist(scale(as.matrix(wsp_tri))), freq_bi_tri,
    comp_bi_tri, sign_bi_tri)

colnames(rect_var) <- c("wrbl_sp", "frequency", "components", "sign")

mfcc_wrbl_mod <- MRM2(formula = wrbl_sp ~ frequency + components +
    sign, nperm = 10000, mat = rect_var)

saveRDS(mfcc_wrbl_mod, "./data/processed/mfcc_warbler_distance_random_start.RDS")

(mfcc_wrbl_mod <- readRDS("./data/processed/mfcc_warbler_distance_random_start.RDS"))

$coef
                wrbl_sp   pval
Int        -0.982621684 0.0001
frequency   1.084976216 0.0001
components  0.077849463 0.0007
sign       -0.003561353 0.6436

$r.squared
       R2      pval 
0.1378832 0.0001000 

$F.test
        F    F.pval 
2295.9832    0.0001 

6.5 warbleR’s fourier cross-correlation

est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")

xc <- cross_correlation(est_schr)

saveRDS(xc, "./data/processed/fourier_cross_correlation_random_start_schroeder.RDS")

xc <- readRDS("./data/processed/fourier_cross_correlation_random_start_schroeder.RDS")

freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(xc),
    "_"), "[[", 2), start = 4, 6)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(xc),
    "_"), "[[", 3), start = 7, 8)))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
    rownames(xc)), "pos", "neg")))

rect_var <- cbind(as.dist(scale(1 - xc)), freq_bi_tri, comp_bi_tri,
    sign_bi_tri)

colnames(rect_var) <- c("fourier_xc", "frequency", "components", "sign")

xc_mod <- MRM2(formula = fourier_xc ~ frequency + components + sign,
    nperm = 10000, mat = rect_var)

saveRDS(xc_mod, "./data/processed/matrix_correlation_fourier_cross_correlation_random_start.RDS")

(xc_mod <- readRDS("./data/processed/matrix_correlation_fourier_cross_correlation_random_start.RDS"))

$coef
            fourier_xc   pval
Int        -1.78485496 0.0001
frequency   1.98172165 0.0001
components  0.14076790 0.0001
sign       -0.02157745 0.0081

$r.squared
       R2      pval 
0.5070224 0.0001000 

$F.test
         F     F.pval 
14764.6553     0.0001 

6.6 Raven’s fourier cross-correlation

xc_rav <- Rraven::imp_corr_mat(file = "random_start_schroeders_raven_fourier_cross_correlation.txt",
    path = "./data/processed")[[1]]

freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(xc_rav),
    "_"), "[[", 2), start = 4, 6)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(xc_rav),
    "_"), "[[", 3), start = 7, 8)))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
    rownames(xc_rav)), "pos", "neg")))

rect_var <- cbind(as.dist(scale(1 - xc_rav)), freq_bi_tri, comp_bi_tri,
    sign_bi_tri)

colnames(rect_var) <- c("fourier_xc_rav", "frequency", "components",
    "sign")

xc_rav_mod <- MRM2(formula = fourier_xc_rav ~ frequency + components +
    sign, nperm = 10000, mat = rect_var)

saveRDS(xc_rav_mod, "./data/processed/matrix_correlation_raven_fourier_cross_correlation_random_start.RDS")

(xc_rav_mod <- readRDS("./data/processed/matrix_correlation_raven_fourier_cross_correlation_random_start.RDS"))

$coef
           fourier_xc_rav   pval
Int            -1.7370728 0.0001
frequency       1.7351569 0.0001
components      0.2531801 0.0001
sign            0.0246393 0.0072

$r.squared
       R2      pval 
0.3757609 0.0001000 

$F.test
        F    F.pval 
8700.3890    0.0001 

6.7 Raven’s waveform correlation

wav_cr_rav <- Rraven::imp_corr_mat(file = "random_start_schroeders_raven_waveform_correlation.txt",
    path = "./data/processed")[[1]]

freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(wav_cr_rav),
    "_"), "[[", 2), start = 4, 6)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(wav_cr_rav),
    "_"), "[[", 3), start = 7, 8)))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
    rownames(wav_cr_rav)), "pos", "neg")))

rect_var <- cbind(as.dist(scale(1 - wav_cr_rav)), freq_bi_tri, comp_bi_tri,
    sign_bi_tri)

colnames(rect_var) <- c("fourier_wav_cr_rav", "frequency", "components",
    "sign")

wav_cr_rav_mod <- MRM2(formula = fourier_wav_cr_rav ~ frequency +
    components + sign, nperm = 10000, mat = rect_var)

saveRDS(wav_cr_rav_mod, "./data/processed/matrix_correlation_raven_waveform_correlation_random_start.RDS")

(wav_cr_rav_mod <- readRDS("./data/processed/matrix_correlation_raven_waveform_correlation_random_start.RDS"))

$coef
           fourier_wav_cr_rav   pval
Int               -1.70508855 0.0001
frequency          1.88607629 0.0001
components         0.04258164 0.0061
sign               0.11819114 0.0001

$r.squared
       R2      pval 
0.4953648 0.0001000 

$F.test
         F     F.pval 
14188.1430     0.0001 

6.8 warbleR’s mel-frequency cross-correlation

est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")

xc <- cross_correlation(est_schr, type = "mfcc")

saveRDS(xc, "./data/processed/mel_frequency_cross_correlation_schroeder_random_start.RDS")

xc <- readRDS("./data/processed/mel_frequency_cross_correlation_schroeder_random_start.RDS")

freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(xc),
    "_"), "[[", 2), start = 4, 6)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(xc),
    "_"), "[[", 3), start = 7, 8)))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
    rownames(xc)), "pos", "neg")))

rect_var <- cbind(as.dist(scale(1 - xc)), freq_bi_tri, comp_bi_tri,
    sign_bi_tri)


colnames(rect_var) <- c("mel_xc", "frequency", "components", "sign")


mxc_mod <- MRM2(formula = mel_xc ~ frequency + components + sign,
    nperm = 10000, mat = rect_var)

saveRDS(mxc_mod, "./data/processed/matrix_correlation_mel_cross_correlation_random_start.RDS")

(mxc_mod <- readRDS("./data/processed/matrix_correlation_mel_cross_correlation_random_start.RDS"))

$coef
                  mel_xc   pval
Int        -1.2613675324 0.0001
frequency   1.2552094577 0.0001
components  0.2486046157 0.0001
sign        0.0008366943 0.9119

$r.squared
       R2      pval 
0.1794016 0.0001000 

$F.test
        F    F.pval 
3138.4764    0.0001 

6.9 DTW on warbleR’s dominant frequency contours

source("~/Dropbox/R_package_testing/warbleR/R/freq_DTW.R")
source("~/Dropbox/R_package_testing/warbleR/R/internal_functions.R")

est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")

dtw_dists <- freq_DTW(est_schr, type = "dominant", img = FALSE)

saveRDS(dtw_dists, "./data/processed/dtw_distance_dominant_frequency_schroeder_random_start.RDS")

dtw_dists <- readRDS("./data/processed/dtw_distance_dominant_frequency_schroeder_random_start.RDS")

freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(dtw_dists),
    "_"), "[[", 2), start = 4, 6)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(dtw_dists),
    "_"), "[[", 3), start = 7, 8)))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
    rownames(dtw_dists)), "pos", "neg")))

rect_var <- cbind(as.dist(scale(dtw_dists)), freq_bi_tri, comp_bi_tri,
    sign_bi_tri)

colnames(rect_var) <- c("dtw_dists", "frequency", "components", "sign")

dtw_dom_mod <- MRM2(formula = dtw_dists ~ frequency + components +
    sign, nperm = 10000, mat = rect_var)

saveRDS(dtw_dom_mod, "./data/processed/matrix_correlation_dtw_distances_dominant_frequency_contours_random_start.RDS")

(dtw_dom_mod <- readRDS("./data/processed/matrix_correlation_dtw_distances_dominant_frequency_contours_random_start.RDS"))

$coef
              dtw_dists   pval
Int        -0.456171403 0.0001
frequency   0.526558411 0.0001
components  0.157208372 0.0001
sign       -0.004413661 0.5422

$r.squared
        R2       pval 
0.02564736 0.00010000 

$F.test
       F   F.pval 
377.8765   0.0001 

6.10 DTW on Raven’s peak frequency contours

source("~/Dropbox/R_package_testing/warbleR/R/freq_DTW.R")
source("~/Dropbox/R_package_testing/warbleR/R/internal_functions.R")

rav_dat <- imp_raven(path = "./data/processed", files = "random_start_schroeder_master_raven_measurements.txt",
    warbler.format = TRUE, all.data = TRUE)

rav_dat <- rav_dat[grep(pattern = "marker", rav_dat$f0, invert = TRUE),
    ]

# rav_dat <- rav_dat[, grep('contour|PFC Slope', ignore.case =
# TRUE, x = names(rav_dat), invert = TRUE)]

rav_freq_cntr <- extract_ts(rav_dat, ts.column = "Peak Freq Contour (Hz)",
    equal.length = TRUE)


est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")

rav_dtw_dists <- freq_DTW(est_schr, ts.df = rav_freq_cntr, type = "dominant",
    img = FALSE)

rownames(rav_dtw_dists) <- rav_dat$orig.sound.file

saveRDS(rav_dtw_dists, "./data/processed/dtw_distance_raven_dominant_frequency_random_start_schroeder.RDS")

rav_dtw_dists <- readRDS("./data/processed/dtw_distance_raven_dominant_frequency_random_start_schroeder.RDS")

freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(rav_dtw_dists),
    "_"), "[[", 2), start = 4, 6)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(rav_dtw_dists),
    "_"), "[[", 3), start = 7, 8)))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
    rownames(rav_dtw_dists)), "pos", "neg")))

rect_var <- cbind(as.dist(scale(rav_dtw_dists)), freq_bi_tri, comp_bi_tri,
    sign_bi_tri)

colnames(rect_var) <- c("rav_dtw_dists", "frequency", "components",
    "sign")

source("~/Dropbox/Projects/geographic_call_variation_yellow-naped_amazon/scripts/MRM2.R")

dtw_rav_peak_mod <- MRM2(formula = rav_dtw_dists ~ frequency + components +
    sign, nperm = 10000, mat = rect_var)

saveRDS(dtw_rav_peak_mod, "./data/processed/matrix_correlation_dtw_distances_raven_dominant_frequency_contours_random_start.RDS")

(dtw_rav_peak_mod <- readRDS("./data/processed/matrix_correlation_dtw_distances_raven_dominant_frequency_contours_random_start.RDS"))

$coef
           rav_dtw_dists   pval
Int        -0.8465155615 0.0001
frequency   1.1803304963 0.0001
components  0.0468889213 0.0323
sign        0.0009260393 0.9175

$r.squared
       R2      pval 
0.1363664 0.0001000 

$F.test
        F    F.pval 
2266.7382    0.0001 

6.11 SAP similarity

• SAP: sound analysis pro
• This analysis was done on 200 ms clips as longer clips did not run in SAP

simil_sap <- readxl::read_xlsx("./data/raw/200ms_SAPsim.xlsx")

# fix weird similarities
simil_sap$`%Similarity`[simil_sap$Accuracy == -1] <- mean(simil_sap$`%Similarity`,
    na.rm = TRUE)

# number of pairwise comparisons posibles
clips <- unique(c(simil_sap$`Sound 1`, simil_sap$`Sound 2`))

ncol(poss_combs <- combn(clips, 2))

poss_combs <- t(poss_combs)

poss_dyad <- sapply(seq_len(nrow(poss_combs)), function(x) {
    paste(sort(c(poss_combs[x, 1], poss_combs[x, 2])), collapse = "|")
})


simil_sap$dyad <- sapply(seq_len(nrow(simil_sap)), function(x) {
    paste(sort(c(simil_sap$`Sound 1`[x], simil_sap$`Sound 2`[x])),
        collapse = "|")
})

# remove duplicated dyads
simil_sap <- simil_sap[!duplicated(simil_sap$dyad), ]

# remove self comparisons
simil_sap <- simil_sap[!simil_sap$`Sound 1` == simil_sap$`Sound 2`,
    ]

nrow(simil_sap)

# missing dyads
missing <- poss_dyad[!poss_dyad %in% simil_sap$dyad]

# split missing by |
missing <- unique(c(sapply(strsplit(missing, "|", fixed = TRUE), "[[",
    1), sapply(strsplit(missing, "|", fixed = TRUE), "[[", 2)))

# remove missing clips
simil_sap <- simil_sap[!(simil_sap$`Sound 1` %in% missing | simil_sap$`Sound 2` %in%
    missing), ]

# get symmetric triangular matrix
dist_sap <- rectangular_to_triangular(as.data.frame(simil_sap))

freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(dist_sap),
    "_"), "[[", 2), start = 4, 6)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(dist_sap),
    "_"), "[[", 3), start = 7, 8)))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
    rownames(dist_sap)), "pos", "neg")))

rect_var <- cbind(as.dist(scale(1 - dist_sap)), freq_bi_tri, comp_bi_tri,
    sign_bi_tri)

colnames(rect_var) <- c("dist_sap", "frequency", "components", "sign")

rect_var[is.na(rect_var[, 1]), 1] <- mean(rect_var[, 1], na.rm = TRUE)

dist_sap_mod <- MRM2(formula = dist_sap ~ frequency + components +
    sign, nperm = 10000, mat = rect_var)

saveRDS(dist_sap_mod, "./data/processed/matrix_correlation_sound_analysis_pro_200_ms_clips.RDS")

(dist_sap_mod <- readRDS("./data/processed/matrix_correlation_sound_analysis_pro_200_ms_clips.RDS"))

$coef
                dist_sap   pval
Int        -1.0474374027 0.0001
frequency   1.3171980171 0.0001
components  0.0285901529 0.1798
sign        0.0006323237 0.9344

$r.squared
      R2     pval 
0.178686 0.000100 

$F.test
        F    F.pval 
3123.2347    0.0001 

$n
[1] 294

6.12 SAP accuracy

• SAP: sound analysis pro
• This analysis was done on 200 ms clips as longer clips did not run in SAP

simil_sap <- readxl::read_xlsx("./data/raw/200ms_SAPsim.xlsx")

# fix weird similarities
simil_sap$Accuracy[simil_sap$Accuracy == -1] <- mean(simil_sap$Accuracy[simil_sap$Accuracy !=
    -1], na.rm = TRUE)

# number of pairwise comparisons posibles
clips <- unique(c(simil_sap$`Sound 1`, simil_sap$`Sound 2`))

ncol(poss_combs <- combn(clips, 2))

poss_combs <- t(poss_combs)

poss_dyad <- sapply(seq_len(nrow(poss_combs)), function(x) {
    paste(sort(c(poss_combs[x, 1], poss_combs[x, 2])), collapse = "|")
})


simil_sap$dyad <- sapply(seq_len(nrow(simil_sap)), function(x) {
    paste(sort(c(simil_sap$`Sound 1`[x], simil_sap$`Sound 2`[x])),
        collapse = "|")
})

# remove duplicated dyads
simil_sap <- simil_sap[!duplicated(simil_sap$dyad), ]

# remove self comparisons
simil_sap <- simil_sap[!simil_sap$`Sound 1` == simil_sap$`Sound 2`,
    ]

nrow(simil_sap)

# missing dyads
missing <- poss_dyad[!poss_dyad %in% simil_sap$dyad]

# split missing by |
missing <- unique(c(sapply(strsplit(missing, "|", fixed = TRUE), "[[",
    1), sapply(strsplit(missing, "|", fixed = TRUE), "[[", 2)))

# remove missing clips
simil_sap <- simil_sap[!(simil_sap$`Sound 1` %in% missing | simil_sap$`Sound 2` %in%
    missing), ]

# get symmetric triangular matrix
dist_sap <- rectangular_to_triangular(as.data.frame(simil_sap[, c("Sound 1",
    "Sound 2", "Accuracy")]))

freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(dist_sap),
    "_"), "[[", 2), start = 4, 6)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(dist_sap),
    "_"), "[[", 3), start = 7, 8)))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
    rownames(dist_sap)), "pos", "neg")))

rect_var <- cbind(as.dist(scale(1 - dist_sap)), freq_bi_tri, comp_bi_tri,
    sign_bi_tri)

colnames(rect_var) <- c("dist_sap", "frequency", "components", "sign")

rect_var[is.na(rect_var[, 1]), 1] <- mean(rect_var[, 1], na.rm = TRUE)

acc_sap_mod <- MRM2(formula = dist_sap ~ frequency + components +
    sign, nperm = 10000, mat = rect_var)

saveRDS(acc_sap_mod, "./data/processed/matrix_correlation_sound_analysis_pro_accuracy_200_ms_clips.RDS")

(acc_sap_mod <- readRDS("./data/processed/matrix_correlation_sound_analysis_pro_accuracy_200_ms_clips.RDS"))

$coef
               dist_sap   pval
Int        -0.791591414 0.0001
frequency   0.801166273 0.0001
components  0.057881918 0.0028
sign        0.003937408 0.5814

$r.squared
       R2      pval 
0.1033568 0.0001000 

$F.test
        F    F.pval 
1654.7898    0.0001 

$n
[1] 294

6.13 SAP harmonic goodness

• SAP: sound analysis pro
• This analysis was done on 200 ms clips as longer clips did not run in SAP

simil_sap <- readxl::read_xlsx("./data/raw/200ms_SAPsim.xlsx")


# number of pairwise comparisons posibles
clips <- unique(c(simil_sap$`Sound 1`, simil_sap$`Sound 2`))

ncol(poss_combs <- combn(clips, 2))

poss_combs <- t(poss_combs)

poss_dyad <- sapply(seq_len(nrow(poss_combs)), function(x) {
    paste(sort(c(poss_combs[x, 1], poss_combs[x, 2])), collapse = "|")
})


simil_sap$dyad <- sapply(seq_len(nrow(simil_sap)), function(x) {
    paste(sort(c(simil_sap$`Sound 1`[x], simil_sap$`Sound 2`[x])),
        collapse = "|")
})

# remove duplicated dyads
simil_sap <- simil_sap[!duplicated(simil_sap$dyad), ]

# remove self comparisons
simil_sap <- simil_sap[!simil_sap$`Sound 1` == simil_sap$`Sound 2`,
    ]

nrow(simil_sap)

# missing dyads
missing <- poss_dyad[!poss_dyad %in% simil_sap$dyad]

# split missing by |
missing <- unique(c(sapply(strsplit(missing, "|", fixed = TRUE), "[[",
    1), sapply(strsplit(missing, "|", fixed = TRUE), "[[", 2)))

# remove missing clips
simil_sap <- simil_sap[!(simil_sap$`Sound 1` %in% missing | simil_sap$`Sound 2` %in%
    missing), ]

# get symmetric triangular matrix
dist_sap <- rectangular_to_triangular(as.data.frame(simil_sap[, c("Sound 1",
    "Sound 2", "Goodness diff")]))

freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(dist_sap),
    "_"), "[[", 2), start = 4, 6)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(dist_sap),
    "_"), "[[", 3), start = 7, 8)))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
    rownames(dist_sap)), "pos", "neg")))

rect_var <- cbind(as.dist(scale(1 - dist_sap)), freq_bi_tri, comp_bi_tri,
    sign_bi_tri)

colnames(rect_var) <- c("dist_sap", "frequency", "components", "sign")

rect_var[is.na(rect_var[, 1]), 1] <- mean(rect_var[, 1], na.rm = TRUE)

good_sap_mod <- MRM2(formula = dist_sap ~ frequency + components +
    sign, nperm = 10000, mat = rect_var)

saveRDS(good_sap_mod, "./data/processed/matrix_correlation_sound_analysis_pro_goodness_200_ms_clips.RDS")

(good_sap_mod <- readRDS("./data/processed/matrix_correlation_sound_analysis_pro_goodness_200_ms_clips.RDS"))

$coef
               dist_sap   pval
Int         0.711874155 0.0001
frequency  -1.182727718 0.0001
components  0.009919279 0.5938
sign        0.000160463 0.9854

$r.squared
       R2      pval 
0.1197922 0.0001000 

$F.test
        F    F.pval 
1953.7405    0.0001 

$n
[1] 294

6.14 SAP pitch difference

• SAP: sound analysis pro
• This analysis was done on 200 ms clips as longer clips did not run in SAP

simil_sap <- readxl::read_xlsx("./data/raw/200ms_SAPsim.xlsx")

# number of pairwise comparisons posibles
clips <- unique(c(simil_sap$`Sound 1`, simil_sap$`Sound 2`))

ncol(poss_combs <- combn(clips, 2))

poss_combs <- t(poss_combs)

poss_dyad <- sapply(seq_len(nrow(poss_combs)), function(x) {
    paste(sort(c(poss_combs[x, 1], poss_combs[x, 2])), collapse = "|")
})


simil_sap$dyad <- sapply(seq_len(nrow(simil_sap)), function(x) {
    paste(sort(c(simil_sap$`Sound 1`[x], simil_sap$`Sound 2`[x])),
        collapse = "|")
})

# remove duplicated dyads
simil_sap <- simil_sap[!duplicated(simil_sap$dyad), ]

# remove self comparisons
simil_sap <- simil_sap[!simil_sap$`Sound 1` == simil_sap$`Sound 2`,
    ]

nrow(simil_sap)

# missing dyads
missing <- poss_dyad[!poss_dyad %in% simil_sap$dyad]

# split missing by |
missing <- unique(c(sapply(strsplit(missing, "|", fixed = TRUE), "[[",
    1), sapply(strsplit(missing, "|", fixed = TRUE), "[[", 2)))

# remove missing clips
simil_sap <- simil_sap[!(simil_sap$`Sound 1` %in% missing | simil_sap$`Sound 2` %in%
    missing), ]

# get symmetric triangular matrix
dist_sap <- rectangular_to_triangular(as.data.frame(simil_sap[, c("Sound 1",
    "Sound 2", "Pitch diff")]))

freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(dist_sap),
    "_"), "[[", 2), start = 4, 6)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(dist_sap),
    "_"), "[[", 3), start = 7, 8)))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
    rownames(dist_sap)), "pos", "neg")))

rect_var <- cbind(as.dist(scale(1 - dist_sap)), freq_bi_tri, comp_bi_tri,
    sign_bi_tri)

colnames(rect_var) <- c("dist_sap", "frequency", "components", "sign")

rect_var[is.na(rect_var[, 1]), 1] <- mean(rect_var[, 1], na.rm = TRUE)

pitch_diff_mod <- MRM2(formula = dist_sap ~ frequency + components +
    sign, nperm = 10000, mat = rect_var)

saveRDS(pitch_diff_mod, "./data/processed/matrix_correlation_sound_analysis_pro_pitch_diff_200_ms_clips.RDS")

(pitch_diff_mod <- readRDS("./data/processed/matrix_correlation_sound_analysis_pro_pitch_diff_200_ms_clips.RDS"))

$coef
              dist_sap   pval
Int         0.26814278 0.0001
frequency  -0.12204666 0.0001
components -0.10396671 0.0001
sign       -0.01128282 0.0745

$r.squared
         R2        pval 
0.003378295 0.000100000 

$F.test
       F   F.pval 
48.66207  0.00010 

$n
[1] 294

6.15 SAP amplitude modulation

• SAP: sound analysis pro
• This analysis was done on 200 ms clips as longer clips did not run in SAP

simil_sap <- readxl::read_xlsx("./data/raw/200ms_SAPsim.xlsx")

# number of pairwise comparisons posibles
clips <- unique(c(simil_sap$`Sound 1`, simil_sap$`Sound 2`))

ncol(poss_combs <- combn(clips, 2))

poss_combs <- t(poss_combs)

poss_dyad <- sapply(seq_len(nrow(poss_combs)), function(x) {
    paste(sort(c(poss_combs[x, 1], poss_combs[x, 2])), collapse = "|")
})


simil_sap$dyad <- sapply(seq_len(nrow(simil_sap)), function(x) {
    paste(sort(c(simil_sap$`Sound 1`[x], simil_sap$`Sound 2`[x])),
        collapse = "|")
})

# remove duplicated dyads
simil_sap <- simil_sap[!duplicated(simil_sap$dyad), ]

# remove self comparisons
simil_sap <- simil_sap[!simil_sap$`Sound 1` == simil_sap$`Sound 2`,
    ]

nrow(simil_sap)

# missing dyads
missing <- poss_dyad[!poss_dyad %in% simil_sap$dyad]

# split missing by |
missing <- unique(c(sapply(strsplit(missing, "|", fixed = TRUE), "[[",
    1), sapply(strsplit(missing, "|", fixed = TRUE), "[[", 2)))

# remove missing clips
simil_sap <- simil_sap[!(simil_sap$`Sound 1` %in% missing | simil_sap$`Sound 2` %in%
    missing), ]

# get symmetric triangular matrix
dist_sap <- rectangular_to_triangular(as.data.frame(simil_sap[, c("Sound 1",
    "Sound 2", "AM diff")]))

freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(dist_sap),
    "_"), "[[", 2), start = 4, 6)))

comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(dist_sap),
    "_"), "[[", 3), start = 7, 8)))

sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
    rownames(dist_sap)), "pos", "neg")))

rect_var <- cbind(as.dist(scale(1 - dist_sap)), freq_bi_tri, comp_bi_tri,
    sign_bi_tri)

colnames(rect_var) <- c("dist_sap", "frequency", "components", "sign")

rect_var[is.na(rect_var[, 1]), 1] <- mean(rect_var[, 1], na.rm = TRUE)

am_sap_mod <- MRM2(formula = dist_sap ~ frequency + components + sign,
    nperm = 10000, mat = rect_var)

saveRDS(am_sap_mod, "./data/processed/matrix_correlation_sound_analysis_pro_amplitude_modulation_200_ms_clips.RDS")

(am_sap_mod <- readRDS("./data/processed/matrix_correlation_sound_analysis_pro_amplitude_modulation_200_ms_clips.RDS"))

$coef
                dist_sap   pval
Int        -0.0239282366 0.2749
frequency   0.0160426132 0.0124
components -0.0075928468 0.4568
sign        0.0008302052 0.8101

$r.squared
          R2         pval 
3.116584e-05 6.990000e-02 

$F.test
        F    F.pval 
0.4474203 0.0699000 

$n
[1] 294

6.16 Combined results

mods <- list(dtw_wv_mod = dtw_wv_mod, wrbl_mod = wrbl_mod, mfcc_wrbl_mod = mfcc_wrbl_mod,
    rav_mod = rav_mod, xc_mod = xc_mod, mxc_mod = mxc_mod, dtw_dom_mod = dtw_dom_mod,
    dtw_rav_peak_mod = dtw_rav_peak_mod, wav_cr_rav_mod = wav_cr_rav_mod,
    xc_rav_mod = xc_rav_mod, dist_sap_mod = dist_sap_mod, acc_sap_mod = acc_sap_mod,
    good_sap_mod = good_sap_mod, pitch_diff_mod = pitch_diff_mod,
    am_sap_mod = am_sap_mod)

names(mods) <- c("DTW waveform period", "warbleR's spectro features",
    "warbleR's MFCC features", "Raven's spectro features", "warbleR's Fourier cross-correlation",
    "Mel cross-correlation", "DTW warbleR's dominant freq", "DTW Raven's peak freq",
    "waveform correlation", "Raven's Fourier cross-correlation", "SAP similarity",
    "SAP accuracy", "SAP goodness diff", "SAP pitch diff", "SAP am")

estimates <- do.call(rbind, lapply(seq_along(mods), function(x) {
    Y <- data.frame(mods[[x]]$coef[-1, ])
    Y$rel_coef <- Y[, 1]/max(Y[, 1])
    Y$mod <- names(mods)[x]
    Y$predictor <- rownames(Y)
    names(Y) <- c("coef", "p", "rel_coef", "model", "predictor")
    return(Y)
}))

estimates$rel_coef <- ifelse(estimates$p < 0.05, estimates$rel_coef,
    0)
estimates$signif <- ifelse(estimates$p < 0.05, "p < 0.05", "p >= 0.05")

estimates$model <- factor(estimates$model, levels = rev(c("SAP pitch diff",
    "SAP am", "SAP goodness diff", "SAP accuracy", "SAP similarity",
    "Raven's spectro features", "Raven's Fourier cross-correlation",
    "waveform correlation", "DTW Raven's peak freq", "warbleR's spectro features",
    "warbleR's MFCC features", "warbleR's Fourier cross-correlation",
    "Mel cross-correlation", "DTW warbleR's dominant freq", "DTW waveform period")))

estimates$predictor <- gsub("sign", "phase", estimates$predictor)

ggplot(estimates, aes(x = predictor, y = model, fill = coef)) + geom_tile() +
    coord_equal() + scale_fill_gradient2(low = viridis(10)[3], high = viridis(10)[7],
    guide = "none") + geom_text(aes(label = round(coef, 3), color = signif),
    size = 2.5) + scale_color_manual(values = c("black", "gray")) +
    labs(x = "", y = "", color = "P value") + theme_classic() + theme(axis.text.x = element_text(color = "black",
    size = 11, angle = 30, vjust = 0.8, hjust = 0.8), axis.text.y = element_text(color = c(rep(viridis(10)[4],
    6), rep(viridis(10)[6], 4))))

Warning: Vectorized input to `element_text()` is not officially supported.
ℹ Results may be unexpected or may change in future versions of ggplot2.

Effect sizes per model and predictor (color intensity shows effect size magnitude relative to the highest effect size within the model).

6.17 Synthesize Shroeders with 2 repetitions

Each one has two repetions follow and are surrounded by 2 100 kHz sine waves that alltogether (sine wave + schroeder + sine wave) add up to a 0.5 s signal

est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")

mean_schroeders <- warbleR:::pblapply_wrblr_int(cl = 20, seq_len(nrow(est_schr)),
    function(x) {
        wave <- read_wave(est_schr, index = x)
        seg <- try_na(mean_segment(wave, plot = FALSE, mean = FALSE,
            type = "ac", thinning = 1))

        return(seg)
    })


names(mean_schroeders) <- paste(est_schr$f0, est_schr$label, sep = "-")

sapply(mean_schroeders, nrow)

mean_schroeders <- mean_schroeders[!sapply(mean_schroeders, function(x) is.na(x[[1]][1]))]


mean_schroeders_list <- lapply(seq_len(length(mean_schroeders)), function(x) {
    data.frame(schroeder = names(mean_schroeders)[x], time = seq(0,
        1, length.out = nrow(mean_schroeders[[x]])), mean.amp = rowMeans(mean_schroeders[[x]]))
})

names(mean_schroeders_list) <- names(mean_schroeders)


file_names <- paste0(f_pad_zero(x = 1:length(mean_schroeders_list),
    width = 3), "_f0-", substr(names(mean_schroeders_list), 0, 3),
    "_ncomp-", sapply(strsplit(gsub("_n|_p", "", names(mean_schroeders_list)),
        "-"), "[[", 2), "_", ifelse(grepl("_p", names(mean_schroeders_list)),
        "pos", "neg"), ".wav")


out <- warbleR:::pblapply_wrblr_int(seq_along(mean_schroeders_list),
    function(x) {

        mean_shc <- mean_schroeders_list[[x]]$mean.amp

        # duplicate it
        mean_shc <- c(mean_shc, mean_shc)

        start_sine <- sine(freq = 100, duration = 44100)@left

        cut_start <- which.min(abs(start_sine[22051:44100] - mean_shc[1])) +
            22050

        start_sine <- start_sine[(cut_start - ((44100/4) - length(mean_shc) -
            1)):cut_start]

        # start_sine <- sine(freq = 100, duration = (44100 / 4)
        # - length(mean_shc))@left

        length(start_sine)

        end_sine <- sine(freq = 100, duration = 44100/2)@left

        cut_end <- which.min(abs(end_sine[1:1000] - mean_shc[length(mean_shc)]))

        end_sine <- end_sine[cut_end:((44100/4) + cut_end - 1)]

        out <- c(start_sine, mean_shc, end_sine)

        wave_obj <- tuneR::normalize(tuneR::Wave(out, samp.rate = 44100,
            bit = 16), unit = "16")

        writeWave(object = wave_obj, filename = file.path("./data/processed/2_repeats_schroeders",
            file_names[x]), extensible = FALSE)

    })


est_rep_schr <- readRDS("./data/processed/extended_selection_table_schroeders.RDS")

# fix start
est_schr <- warbleR::selection_table(whole.recs = TRUE, path = "./data/processed/2_repeats_schroeders",
    extended = TRUE, confirm.extended = FALSE)

est_schr$bottom.freq <- est_rep_schr$bottom.freq

est_schr$top.freq <- est_rep_schr$top.freq

est_schr$f0 <- gsub("f0-", "", sapply(strsplit(est_schr$sound.files,
    split = "_"), "[[", 2))

est_schr$label <- gsub("os.wav_1|eg.wav_1", "", sapply(strsplit(est_schr$sound.files,
    split = "ncomp-"), "[[", 2))

master <- master_sound_file(X = est_schr, file.name = "schroeder_master_2_repeats",
    dest.path = "./data/processed/", gap.duration = 0.8, cex = 14)

write.csv(master, "./data/processed/master_annotations_schroeders_2_repeats.csv",
    row.names = FALSE)

master$f0 <- c("start_marker", paste("f0", est_schr$f0), "")

master$label <- c("marker", est_schr$label, "marker")

Rraven::exp_raven(master, path = "./data/processed/", sound.file.path = "./data/processed/",
    file.name = "master_annotations_schroeders_2_repeats_raven")

warbleR::full_spectrograms(master, path = "./data/processed/", sxrow = duration(readWave("./data/processed/schroeder_master_2_repeats.wav"))/12,
    rows = 12, horizontal = TRUE, dest.path = "./data/processed",
    collevels = seq(-150, 0, 5), labels = "label", song = "f0", wl = 2000,
    suffix = "2_repeats")

master sound file

 

7 Takeaways

• Amplitude autocorrelation works better at getting periodicity
• Dynamic time warping dissimilarity on mean periods, Raven’s waveform correlation and (to some extent) Fourier cross-correlation did capture phase differences in Schroeders
• Dynamic time warping dissimilarity on mean periods remove differences in frequency when forcing all periods to have the same number of values (i.e. same length)

 

---

 

Session information

R version 4.4.2 (2024-10-31)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Costa_Rica
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] animation_2.7        numform_0.7.0        ecodist_2.1.3       
 [4] PhenotypeSpace_0.1.0 ggplot2_3.5.1        baRulho_2.1.4       
 [7] ohun_1.0.2           Rraven_1.0.14        viridis_0.6.5       
[10] viridisLite_0.4.2    warbleR_1.1.34       NatureSounds_1.0.5  
[13] tuneR_1.4.7          seewave_2.2.3        formatR_1.14        
[16] knitr_1.49           kableExtra_1.4.0    

loaded via a namespace (and not attached):
 [1] DBI_1.2.3              bitops_1.0-9           deldir_2.0-4          
 [4] pbapply_1.7-2          gridExtra_2.3          remotes_2.5.0         
 [7] permute_0.9-7          testthat_3.2.3         rlang_1.1.5           
[10] magrittr_2.0.3         e1071_1.7-16           compiler_4.4.2        
[13] spatstat.geom_3.3-2    mgcv_1.9-1             png_0.1-8             
[16] systemfonts_1.1.0      vctrs_0.6.5            Sim.DiffProc_4.9      
[19] stringr_1.5.1          pkgconfig_2.0.3        crayon_1.5.3          
[22] fastmap_1.2.0          backports_1.5.0        labeling_0.4.3        
[25] rmarkdown_2.28         xfun_0.50              jsonlite_1.8.9        
[28] goftest_1.2-3          spatstat.utils_3.1-0   terra_1.7-78          
[31] Deriv_4.1.6            parallel_4.4.2         cluster_2.1.6         
[34] R6_2.6.1               stringi_1.8.4          spatstat.data_3.1-2   
[37] spatstat.univar_3.0-1  brio_1.1.5             Rcpp_1.0.14           
[40] tensor_1.5             splines_4.4.2          Matrix_1.7-0          
[43] igraph_2.1.4           tidyselect_1.2.1       rstudioapi_0.16.0     
[46] abind_1.4-8            yaml_2.3.10            vegan_2.6-8           
[49] codetools_0.2-20       dtw_1.23-1             spatstat.random_3.3-1 
[52] spatstat.explore_3.3-2 curl_6.2.0             lattice_0.22-6        
[55] tibble_3.2.1           withr_3.0.2            evaluate_1.0.3        
[58] signal_1.8-1           sf_1.0-19              sketchy_1.0.5         
[61] units_0.8-5            proxy_0.4-27           polyclip_1.10-7       
[64] xml2_1.3.6             pillar_1.10.1          packrat_0.9.2         
[67] KernSmooth_2.23-24     checkmate_2.3.2        generics_0.1.3        
[70] sp_2.1-4               RCurl_1.98-1.16        munsell_0.5.1         
[73] scales_1.3.0           class_7.3-22           glue_1.8.0            
[76] tools_4.4.2            xaringanExtra_0.8.0    grid_4.4.2            
[79] colorspace_2.1-1       nlme_3.1-165           raster_3.6-26         
[82] cli_3.6.4              spatstat.sparse_3.1-0  svglite_2.1.3         
[85] dplyr_1.1.4            gtable_0.3.6           fftw_1.0-9            
[88] digest_0.6.37          classInt_0.4-11        farver_2.1.2          
[91] rjson_0.2.23           htmlwidgets_1.6.4      htmltools_0.5.8.1     
[94] lifecycle_1.0.4        httr_1.4.7             MASS_7.3-61