Audio segmentation using Flattened Local Trimmed Range for ecological acoustic space analysis

Modeling the spectrogram

We model the spectrogram S_db(t, f) as a sum of different energetic processes: (1)${\displaystyle S_{db}\left(t,f\right)=b\left(f\right)+\epsilon \left(t,f\right)+\sum _{i=1}^n{R}_i\left(t,f\right){\mathbb{I}}_i\left(t,f\right),}$ where b(f) is a frequency-dependent process that is taken as constant in time, while ϵ(t, f) is a process that is stationary in time and frequency with 0-mean, 0-median, and some scale parameter, and R_i(t, f) is a set of non-zero-mean localized energy processes that are bounded by their support functions ${\mathbb{I}}_i\left(t,f\right)$, for 1 ≤ i ≤ n. An interpretation for these energetic processes is that b(f) corresponds to a frequency-dependent near-constant noise, ϵ(t, f) corresponds to a global noise process with a symmetric distribution and the R_i(t, f) are our acoustic events, of which there are n.

In this model, we assume that the set of localized energy processes has four properties:

• A1

  The localized energy processes are mutually exclusive and are not adjacent. That is, no two localized energy processes share in the same (t, f) coordinate, nor do they have adjacent coordinates. Thus, ∀1 ≤ i, j ≤ n, i ≠ j, 0 ≤ t < τ, 0 ≤ f < η, we have: ${\displaystyle {\mathbb{I}}_i\left(t,f\right){\mathbb{I}}_j\left(t,f\right)=0}$ ${\displaystyle {\mathbb{I}}_i\left(t+1,f\right){\mathbb{I}}_j\left(t,f\right)=0}$ ${\displaystyle {\mathbb{I}}_i\left(t-1,f\right){\mathbb{I}}_j\left(t,f\right)=0}$ ${\displaystyle {\mathbb{I}}_i\left(t,f+1\right){\mathbb{I}}_j\left(t,f\right)=0}$ ${\displaystyle {\mathbb{I}}_i\left(t,f-1\right){\mathbb{I}}_j\left(t,f\right)=0.}$

  This can be done without loss of generality. If two such localized processes existed, we can just consider their union as one.

• A2

  The regions of localized energy processes dominate the energy distribution in the spectrogram on each given band. That is, ∀0 ≤ t₁, t₂ < τ, 0 ≤ f < η, we have: ${\displaystyle ϵ\left(t_1,f\right)+b\left(f\right)\le ϵ\left(t_2,f\right)+b\left(f\right)+\sum _{i=1}^n{\mathbb{I}}_i\left(t_2,f\right)R_i\left(t_2,f\right).}$

• A3

  The proportion of samples within a localized energy processes in a given frequency band, denoted as ρ(f), is less than half the samples in the entire frequency band. That is, ∀0 ≤ f < η, we have: ${\displaystyle \rho \left(f\right)=\frac{1}{\tau }\sum _{t=0}^{\tau }\sum _{i=1}^n{\mathbb{I}}_i\left(t,f\right)<.5.}$

• A4

  Each localized energy process dominates the energy distribution in its surrounding region, when accounting for frequency band-dependent effects. That is, for every (t₁, f₁) point that falls inside a localized energy process (∀1 ≤ i ≤ n, 0 ≤ t₁ < τ, 0 ≤ f₁ < η) where ${\mathbb{I}}_i\left(t_1,f_1\right)=1$), there is a region-dependent time-based radius r_i,1 and a frequency-based radius r_i,2, such that for every other (t₂, f₂) point around the vicinity (∀(t₂, f₂) ∈ [t₁ − r_i,1, t₁ + r_i,1] × [f₁ − r_i,2, f₁ + r_i,2]), we have: ${\displaystyle ϵ\left(t_2,f_2\right)\le ϵ\left(t_1,f_1\right)+\sum _{i=1}^n{\mathbb{I}}_i\left(t_1,f_1\right)R_i\left(t_1,f_1\right).}$

We want to extract the R_i components or, more specifically, their support functions ${\mathbb{I}}_i\left(t,f\right)$, from the spectrogram. If we are able to estimate b(f) reliably for a given spectrogram, we can then compute $\hat{S}\left(t,f\right)=S_{db}\left(t,f\right)-b\left(f\right)$, a spectrogram that is corrected for frequency intensity variations. Once that is done, we compute local statistics to estimate the ${\mathbb{I}}_i\left(t,f\right)$ regions and, thus, segregate the spectrogram into the localized energy processes R_i(t, f) and an ϵ(t, f) background process.

Flattening—estimating b(f)

Other than A2, A3 and A4, we do not hold any assumptions for ϵ(t, f) or R_i(t, f). In particular we do not presume to know their distributions. Thus, it is difficult to formulate a model to compute a Maximum A-Posteriori Estimator of b(f)|S_db(t, f). Even so, the frequency sample means $\mathrm{\mu }\left(f\right)=\frac{1}{\tau }{\sum }_{t=0}^{\tau -1}{S}_{db}\left(t,f\right)$ of a given spectrogram do not give a good estimate on b(f) since they get mixed with the sum of non-zero expectations of any intersecting region: ${\displaystyle \mathrm{\mu }\left(f\right)=\frac{1}{\tau }\sum _{t=0}^{\tau -1}{S}_{db}\left(t,f\right)=\frac{1}{\tau }\sum _{t=0}^{\tau -1}\left(b\left(f\right)+ϵ\left(t,f\right)+\sum _{i=1}^n{R}_i\left(t,f\right){\mathbb{I}}_i\left(t,f\right)\right)=b\left(f\right)+\frac{1}{\tau }\sum _{t=0}^{\tau -1}\sum _{i=1}^n{R}_i\left(t,f\right){\mathbb{I}}_i\left(t,f\right).}$

Since ϵ(t, f) is a stationary 0-mean process, we do not need to worry about it as it will eventually cancel itself out, but the localized energy process regions do not cancel out. Since our goal is to separate these regions from the rest of the spectrogram in a general manner, if an estimate of b(f) is to be useful, it should not depend on the particular values within these regions.

While using the mean does not prove to be useful, we can use the frequency sample medians, along with A2 and A3 to remove any frequency-dependent time-constant bands from the spectrogram. We formalize this with the following theorem:

Thus, assuming A2 and A3, m(f) gives an estimator whose bias is limited by the range of the ϵ process and is completely unaffected by the R_i processes. Furthermore, as ρ(f) approaches 0, m(f) approaches b(f).

We use the term band flattening to refer to the process of subtracting the b(f) component from S_db(t, f). Thus we call $\hat{S}\left(t,f\right)=S_{db}\left(t,f\right)-m\left(f\right)$ the band flattened spectrogram estimate of S_db(t, f). Figure 2B shows the output from this flattening procedure. As can be seen, this procedure removes any frequency-dependent time-constant bands in the spectrogram.

Estimating ${\mathbb{I}}_i\left(t,f\right)$

We can use the band flattened spectrogram $\hat{S}\left(t,f\right)$ to further estimate the ${\mathbb{I}}_i\left(t,f\right)$ regions, since: ${\displaystyle \hat{S}\left(t,f\right)\approx S_{db}\left(t,f\right)-b\left(f\right)=ϵ\left(t,f\right)+\sum _{i=1}^n{R}_i\left(t,f\right){\mathbb{I}}_i\left(t,f\right).}$

We do this by computing the local α-trimmed range ${\mathrm{Ra}}_{\alpha }\left\{\hat{S}\right\}$. That is, given some 0 ≤ α < 50, and some r > 0, for each (t, f) pair, we compute: ${\displaystyle {\mathrm{Ra}}_{\alpha }\left\{\hat{S}\right\}\left(t,f\right)=P_{100-\alpha }\left({\hat{S}}_{{\mathrm{ne}}_r\left(t,f\right)}\right)-P_{\alpha }\left({\hat{S}}_{{\mathrm{ne}}_r\left(t,f\right)}\right),}$ where P_α(⋅) is the α percentile statistic, and ${\hat{S}}_{{\mathrm{ne}}_r\left(t,f\right)}$ is the band flattened spectrogram, with its domain restricted to a square neighborhood of range r (in time and frequency) around the point (t, f).

Assuming A4, the estimator would give small values for neighborhoods without localized energy processes, but would peak around the borders of any such process. This statistic could then be thresholded to compute estimates of these borders and an estimate of the support functions ${\mathbb{I}}_i\left(t,f\right)$. Figure 2C shows the local trimmed range of a flattened spectrogram image. As can be seen, areas with acoustic events have a higher local trimmed range, while empty areas have a lower one.

FLTR validation

We used the FLTR algorithm with a 21 × 21 window and α = 5 to extract the acoustic events and compared them with the manually labeled acoustic events.

For the validation, we used two comparison methods, the first is based on a basic intersection test between the automatically detected and the manually labeled events’ bounds, and the second one is based on an overlap percent. For each manual label and detection event pair, we defined the computed overlap percent as the ratio of the area of their intersection and the area of their union: ${\displaystyle \mathbb{O}\left(L,D\right)=\frac{\mathbb{A}\left(L\cap D\right)}{A\left(L\cup D\right)},}$ where L is a manually labeled event, D is a automatically detected event area, and $\mathbb{A}\left(L\cap D\right)$ and $\mathbb{A}\left(L\cup D\right)$ are the area of the intersection and the union of their respective bounds.

On the first comparison method, for each acoustic event whose bounds intersected the bounds of a manually labeled acoustic event, we registered it as a detected acoustic event. For each detected event whose bounds did not intersect any manually labeled acoustic events, we registered it as detected, but without an acoustic event. On the other hand, the manually labeled acoustic events that did not intersect any detected acoustic event were registered as undetected acoustic events.

The second method followed a similar path as the first, but it requires an overlap percent of at least 25%. For each acoustic event whose overlap percent with a manually labeled acoustic event was greater than or equal to 25%, we registered it as a detected acoustic event. For each detected event that did not have any manually labeled acoustic events with an overlap percent of at least 25%, we registered it as detected, but without an acoustic event. On the other hand, the manually labeled acoustic events that did not have an overlap percent of at least 25% with any detected acoustic event were registered as undetected acoustic events.

These data were used to create a confusion matrix to compute the FLTR algorithm’s sensitivity and positive predictive value for each method. The sensitivity is computed as the ratio of the number of manually labeled acoustic events that were automatically detected (true positives) over the total count of manually labeled acoustic events (true positives and false negatives). ${\displaystyle \text{Sensitivity}=\frac{\text{True Positives}}{\text{True Positives}+\text{False Negatives}}.}$ This measurement reflects the percent of detected acoustic events that were present in the recording.

The positive predictive value is computed as the ratio of the number of manually labeled acoustic events that were automatically detected (true positives) over the total count of detected acoustic events (true positives and false positives). ${\displaystyle \text{Positive Predictive Value}=\frac{\text{True Positives}}{\text{True Positives}+\text{False Positives}}.}$ This measurement reflects the percent of real acoustic events among the set of detected acoustic events.

FLTR application

As with the FLTR validation step, we used the FLTR algorithm with a 21 × 21 window and α = 5 to extract the acoustic events in each of the recording samples in the second dataset, which we then converted into feature vectors.

The variables computed for each extracted acoustic event R_i were:

• Time of Day. Hour in which the recording from this acoustic event was taken.

• Bandwidth. Length of the acoustic event in Hertz. That is, given $F_i=\left\{f|\exists t,{\mathbb{I}}_i\left(t,f\right)=1\right\}$, the bandwidth is defined as: (2)${\displaystyle {\text{bw}}_i=max\left(F_i\right)-min\left(F_i\right).}$

• Duration. Length of the acoustic event in seconds. That is, given $T_i=\left\{t|\exists f,{\mathbb{I}}_i\left(t,f\right)=1\right\}$, the duration is defined as: (3)${\displaystyle {\text{dur}}_i=max\left(T_i\right)-min\left(T_i\right).}$

• Dominant Frequency. The frequency at which the acoustic event attains its maximum power: (4)${\displaystyle {y_{\text{max}}}_i={argmax}_f\left(\underset{t}{max}{\mathbb{I}}_i\left(t,f\right)\hat{S}\left(t,f\right)\right).}$

• Coverage Ratio. The ratio between the area covered by the detected acoustic event and the area detected by the bounds: (5)${\displaystyle {\text{cov}}_i=\frac{\sum _{t,f}{\mathbb{I}}_i\left(t,f\right)}{{\text{dur}}_i{\text{bw}}_i}.}$

Using these features, we generate a log-density plot matrix for all pairwise combinations of the features for each site. The plots in the diagonal are log histograms of the column feature.

We measured the information content of each feature pair by computing their joint entropy H: ${\displaystyle H=\frac{1}{N_1{N}_2}\sum _{i=1}^{N_1}\sum _{j=1}^{N_2}{h}_{i,j}{log}_2{h}_{i,j},}$ where log₂ is the base 2 logarithm, h_i,j is the number of events in the (i, j)th bin of the joint histogram, and N₁ and N₂ are the number of bins of each variable in the histogram. A higher value of H means a higher information content.

We also focused our attention on areas with a high and medium count of detected acoustic events (log of detected events greater than 6 and 3.5, respectively). As an example, we selected an area of interest in the feature space of the Sabana Seca dataset (i.e., a visual cluster in the bw vs. y˙max plot). We sampled 50 detected acoustic events from the area and categorized them manually.

Sampled areas of interest

In the Sabana Seca bw vs. y˙max plot, there are five high event detection count areas (Fig. 7). We focus on the region with bandwidth between 700 Hz and 1,900 Hz, and dominant frequency between 2,650 Hz and 3,550 Hz.

The 50 sampled acoustic events from the area of interest were arranged into six groups (Table 3). The majority of the events (40 out of 50 events) were the second note, “qui,” of the call of the frog Eleutherodactylus coqui. The second largest group is composed of the chirp sound of an Leptodactylus albilabiris (5 events). The third group is two unknown, similar, calls. The fourth group is an event of an Leptodactylus albilabiris chirp with an almost overlapping Eleutherodactylus coqui’s “qui” tone. The last two groups are the call of an unknown insect, and an acoustic event arising from a person speaking on a radio station from an interference picked up by the recorder.