Identify abrupt changes in an audio waveform

Question

I have an audio file which was created by substituting multiple speech segments into an original utterance. I want to find the timestamps where such substitutions happened, so I guess that one thing I can do is looking for the abrupt changes in the audio waveform.

By visualizing the waveform, I have something like below. There was a sudden shift in the middle of the figure, which implies the fusion of the two speech segments. librosa version is 0.9.2 and the audio sample can be retrieved here.

fig, (ax, ax2) = plt.subplots(nrows=2, sharex=True)
ax.set(xlim=[1.5, 1.6], title='Sample view', ylim=[-1, 1])
y, sr = librosa.load("example.wav", sr=None,duration=10)
y_harm, y_perc = librosa.effects.hpss(y)
librosa.display.waveshow(y, sr=sr, ax=ax, marker='.', label='Full signal')
librosa.display.waveshow(y_harm, sr=sr, alpha=0.5, ax=ax2, label='Harmonic')
librosa.display.waveshow(y_perc, sr=sr, color='r', alpha=0.5, ax=ax2, label='Percussive')

I found some answer here and here which suggested

differentiate the signal numerically and threshold the magnitude of differentiate to identify the abrupt changes

My question is: how to find such abrupt changes and their timestamps programmatically by using any (maybe, python-based) toolkits? Assume that there are only two labels per audio file, each label represents a speaker or an unknown distribution of speech.

Note that the segment merged done by waveform overlapping and the amplitude is normalized. I think it will make it harder to locate such "joints". I am not sure I used to right terms for this. Please suggest the right ones if you I think there are proper terms to denote what I have described so far.

Answer 1

Synchrosqueezed Wavelet Transform is an option. I have developed a complete algorithm for this task, which scores 100% train accuracy and 86% test accuracy with 0.05 sec tolerance, without machine learning, on the provided dataset of 14 samples. Its features can be reused with ML.

I explain my algorithm, suggest improvements, and recommend ML-oriented alternatives.

Overview

"Detection" can be seen as a case of classification, and classification is about balancing discriminativeness and invariance - or simply, maximizing sensitivity to relevant variability. Sparsity and selective activation are excellent properties for this, and SSQ on clean audio excels at this.

The "abrupt change" to detect here is the immediate insertion of a new audio segment, and silencing of previous - a sharp discontinuity. Equivalently, it's the superposition of two audios with zero-padding - and zero-padding is well-studied, by me. It causes predictable distortions in time-frequency, lowering SNR - but for us, this noise is the signal.

Quick inspection: CWT vs STFT

The log scale is favored for many natural signals, especially audio, but check anyway:

CWT richly maps out the variability that's compactified in STFT. Mel-spectrogram is also an option but not explored here.

Inspect regions of interest

We zoom around 1.55 sec and 3.50 sec,

SSQ is sparse when, in a 2D time-frequency slice, complex-valued oscillations share a frequency. Wavelet inner products, faced with a discontinuity, stop having this property - so SSQ fails to remap to a common frequency bin, and is instead spread out smoothly.

Amplify discriminative features / Attenuate invariants

"Noise is signal", so "signal is noise" - so we want to get rid of the signal, which here is transformed into tight curves/lines in 2D.

We need some way to target the ridges (curves); ridge extraction is an option but can be hard to work with. Instead we can exploit the tight concentration of ridges with an otherwise terrible derivative operator - diff, i.e. $x[n] - x[n - 1]$, along the frequency axis, which acts as a highpass filter with extreme time resolution ($h[n] = [1, -1]$). Sharp transitions are high-frequency structures, which will be amplified - and the smooth zero-padding artifacts will be attenuated - all while preserving the original spread of features per high time resolution (not so with e.g. Gaussian highpass filter).

So, diff(Tx, axis=0) where Tx = abs(ssq_cwt(x)):

Now we apply a thresholding upon Tx_diff, to zero points in Tx. For this I use the sparse_mean algorithm I developed - a more robust alternative to mean in that it's unaffected by the pure size of input, so sparse_mean(zero_pad(x)) == sparse_mean(x) while mean(zero_pad(x)) != mean(x). This makes "mean" of Tx of one sine, same as mean for two or more sines, for example.

So, Tx[Tx_diff > sparse_mean(Tx_diff)] = 0, and call it Tx_carved:

Ah-ha! It's red where we want it. But not just there - there's many dots scattered around as thresholding leftovers, that we can't get rid of by re-thresholding. The distinguishing feature here is the size of area that's red, which can be measured by aggregation. Design as follows:

1. There's two simple options - one being 2D filtering, as 2D local sum. But, a key piece of the puzzle is that the impulse response of CWT & STFT is vertically concentrated, and response to impulse resembles response to general spikes and sharp deformations. So, we favor vertical summation over all frequencies.
2. But we notice that the feature is also spread out in time - so we also spread in time, by the expected spread amount, which should approximately be the average size of the time-frequency kernel's temporal envelope.
3. Lastly, we can sum absolute value, or its square (energy). Energy can be shown to accurately reflect "amount of information" (as I discovered), and is also the conserved quantity between time and frequency domains, so we favor it.

All together: (1) sum Tx_carved**2 along frequency; (2) compute running sum by convolving with brickwall - where red vertical lines are our labels (ground truth):

Good enough, I'd say.

Making predictions

The error scores were computed as abs(true_index - pred_index) / len(x). Two kinds of predictions were made:

1. one_peak_pred = peak; peak = argmax(feature_1d), self-explanatory.
2. two_peaks_pred. We zero the interval in feature_1d from peak - pwidth to peak + pwidth, where pwidth is the width of the wavelet computed earlier, and the width of our brickwall. Then we repeat argmax(feature_1d), that's peak2, and take (peak + peak2) / 2.

I'll comment on which one's best later.

The case of one prediction is clear, there's only one argmax - but what about two or more? If we zero peak then take argmax again, the result could be (and almost certainly will be) right next to it. The previous visual has a major lapse, in that predictions were made by first zooming around the labels - the full vector:

Ouch! The top two distinct peaks are both wrong.

Amplify discriminants (2): template matching

Observe Tx_carved again:

We note that the features of interest have a particular shape. Indeed, it's the impulse response of the CWT filterbank, which SSQ completely fails to reassign (as described earlier). So, instead of a naive local energy aggregation, we try convolving with the 2D impulse response - but only once, along time, left to right (not full 2D convolution):

Much better. And the updated zoom:

I call it Subtle Impulse Intensity (SII) vector: "subtle" since its norm in original CWT is much lower than the signal's; "impulse" since its shape is CWT's impulse response; "intensity" since it derives from inner product (similarity) with CWT magnitude (intensity).

Making predictions (2)

By eye, there's clearly two distinct peaks. How to code?

Observe that both peaks have approximately the same width - that's not coincidence, it's the effective temporal width, pwidth, of the 2D impulse response. So we measure this width, and after finding peak, zero the interval around it: SII[peak - pwidth:peak + pwidth] = 0:

The next argmax(SII) will be exactly where we want.

Amplify discriminants (3): second derivative carving

We apply a second derivative, this time upon Tx_carved, then carve it again:

and SII is now

This is something I discovered by accident and haven't fully understood. It appears to clean up the sharp leftovers of the first pass - but further passes do not yield improvements.

Amplify discriminants (4): energy-frequency scaling

Here's the full CWT and SSQ_CWT:

Minimal low frequency content - that's expected. But also, the last two octaves - $2000\ \text{Hz}$ to $4000\ \text{Hz}$ to $8000\ \text{Hz}$ - are less bright. That's not because there's less actual information there - it's the "square power law" of audio: $P \sim 1/f^2$ (or so I've read in a few places... a reference):

Note we must scale via np.logspace() rather than np.arange() per CWT. Applying, we obtain:

Looks ideal.

Except, while I did favor higher frequencies, I didn't quite do it like this - see Addendum.

Attenuate invariants (2): lessen template confounding

I said CWT is better, but it's not better at everything. The algorithm up to this point fails at another example:

It misses the label near 1.5 and favors one around 0.7, so let's inspect both (energy-scaled):

Indeed there's high-frequency signal content that's not sufficiently reassigned by SSQ, so carving ignores it. If we listen to the file, T2_D_00004326.wav, near 0.7, we hear a sharply spoken "tya"; there's also many "sh" later, also high frequency.

SSQ_STFT (energy-scaled):

Much more informative - it's the reverse case compared to our opening example, where STFT richly maps out the frequency interval that's compactified by CWT. Note that STFT can perfectly track LFM and CWT EFM (linear / exponential frequency modulation), and over sufficiently large frequency intervals, it's impossible for one transform to match the performance of the other just by tuning the wavelet/window.

So, we go ahead and carve, this time without the derivative (the ridges aren't as sharp), then invert (via one-integral), and then pass it to our existing algorithm. But here we should be conservative and carve little, as eliminating too much of the signal will also eliminate the resulting impulse. New results with this "STFT pre-filtering":

Still leftovers, but much less. New SII:

Other steps

It's too much for one post, so I summarize:

1. Ignore silence: regions of approximate silence are detected and omitted from (1) generating SII; (2) generating predictions.

   • The minimum contiguous interval to qualify as "ignorable silence" is set to be silence_interval = 0.2, and minimum distance to a prediction as silence_proximity = 0.2, both in seconds. This is based on an estimation as to what comprises the greatest possible interval between "transition" of audios, and is useful for any model, as such an interval must exist, otherwise it's impossible to distinguish a "transition" from simply silence.
   • It was motivated by SII having large peaks due to sharp behavior in silences as audio is just starting; since SII is based on small-norm time-frequency patterns, and silence is small norm by definition, it makes false positives likelier.

2. Duration-based exclusion: I wrote that SII around peak is carved according to pwidth, but that was only originally so. As long as we've not accounted for all invariants, there will be false positives, so it helps to eliminate as large an interval around peak as possible, and I set that to min_audio_interval = 0.5, which is an estimate of the shortest possible audio segment, before or after transition.

3. Multi-stage energy scaling: for various reasons, I found it better to rescale not just the original SSQ_CWT, but also SSQ_STFT and the 2D impulse response, as well as Tx_carved: four scalings in total.

4. Zero-padding up to scoring: I zero pad and don't unpad until it's time to find peak. Zero-padding is the correct prior (assume no audio before or after the segment), and operations between scoring and ssq_cwt are subject to boundary effects (left edge influencing right edge and vice versa).

5. Moving sum of energy of SII: after cross-correlating with the 2D impulse response, I square the result (for reasons described earlier), and convolve with a rectangular window of width pwidth, as a way of generating a vector of local sums of energies over an interval of a relevant size. That's because we expect the most of the impulse feature's energy to be confined within pwidth in SSQ_CWT.

6. Favor one-peak prediction: originally it was two-peak, but after all the refinings, I no longer found any two peaks, and the carving around peak degraded two-peak performance.

All S.I.I. vectors

Performance improvement suggestions

I rewrote the code a bit, now it's 17 sec CPU and 6 sec GPU for a 5 sec file (i7-7700HQ, GTX 1070). This can be made way faster by (1) moving everything to GPU, basically finding torch equivalents to np; (2) batching; (3) intelligently cutting down operations. (1) is easy, (2) tedious, (3) hard.

ELI5

SSQ upon speech makes pictures with sharp curves. Speech transitions are cones in those pictures. Speech is bright but cones are dim, so remove speech by erasing what's sharp. Put a transparent cone over the picture, slide it left to right, and underneath draw a line with a pen; if the picture looks more like that cone, draw higher. Once that's done, the two highest points drawn are our predictions.

Possible followup questions

1. How can SSQ_STFT be used for this task? -- There's more than I described; it managed to produce better discriminants on where my algorithm failed, but also much worse invariants, and I've not had the time to work it out. This approach has the potential to be much faster than the fastest version of my algorithm.

2. How to reliably find mean of sparse data? -- For thresholding or other uses. Outside this answer, I use sparse_mean to find c in log norm: log(1 + x/c).

3. How is 1D cross-correlation with 2D impulse response efficiently computed? -- Actually I already did it.

---

ML recommendations

Success with ML + features requires not only properly configuring the feature extractor with respect to data, but also the ML algorithm with respect to the transform. This is critical, and applicable to all below.

1. Normalize properly: "log norm" almost never means log(x), it's log(1 + x/c), and the right c can make all the difference. Definitely avoid x /= max(abs(x)), favor x /= x.std().
2. Mel-spectrogram: I know little of it besides "take STFT and log-scale it", but that's plenty. OP stated they've tried and failed, I point to my note above.
3. SSQ_CWT: this is a troubling feature to pair with ML but can be done
4. Carved SSQ_CWT: it can either be used as a feature itself, or be inverted as a form of highly nonlinear "de-signaling", then fed to other feature extractors
5. SII: I don't recall seeing any false negatives; failed predictions were only a result of false ranking. Hence it makes for an excellent "attention" vector to feed alongside other features to an NN.
6. Spectrogram: while I spoke against it, recently I found something very convenient that I've not had the time to fully explore
7. ??? This is likely my top recommendation, but for reasons, I can't disclose it at the moment. Will update in the near future if I remember.

Other answers

In short, I doubt any of the suggestions will work on their own, and will require significant further feature engineering. I'll comment on one of the answers right now:

Royi: at first I thought this may work, but not after inspecting the data. Some transitions are difficult to tell even by ear: the speaker, tone, tempo, and kind of words spoken can all be same - so the nonparametric models, which from I understand are anomaly detectors, are unlikely to work. The issue is that the norm of the transition event that constitutes the signal is very low even in time-frequency, so working with raw data will be catastrophic. I know nothing about Kalman but same concern applies.

---

Addendum

Energy vs frequency scaling

Checking the code while I'm writing the post, I notice that my scheme differs tremendously from a direct $1/f^2$ scaling:

1. I scale at four separate stages - three CWT, one STFT. For STFT it's actually exactly $1/f^2$.
2. For CWT, firstly, I manually set scaling limits at 0.1 and 1, rather than according to the actual wavelet frequencies - and that's with respect to the CWT trimmed according to fmax_idx_frac.
3. I squared the scaler that's applied to amplitudes (modulus), not energies (modulus squared).
4. One of the stages involves scaling the impulse response that's convolved (cross-correlated) with, and that was scaled by the square of the square - unintended, coding error.

So, I went ahead and fixed everything - and accuracy dropped (+2 missed predictions), and SII's were much worse. Possible explanations:

1. The "square law" isn't exact nor agnostic to high vs low regimes; see Stokes's law of sound attenuation. For one, ratio of min to max frequency in CWT here was x32,000 - this completely obliterated non-high frequencies, and mind that this is applied thrice. So whatever I did may have simply coincided better with reality.
2. Even if the law were exact, there may be other factors - namely, greater importance of high frequencies per the impulse response being less subject to distortions there. Indeed, before I implemented the STFT high frequency pre-filtering, there were more false positives - so perhaps we also need a low frequency pre-filtering.
3. I've not changed any thresholding while testing. The ones in carving may certainly have required it.
4. (Realized later) I've not actually "fixed everything", I probably should rescale filters according to Littlewood-Paley sum, I've not tried. It doesn't make sense for it to work for STFT but not CWT.

---

Citation

If you use a significant portion of this answer or code, please cite it as

John Muradeli, 2023. Identify abrupt changes in an audio waveform. URL: https://dsp.stackexchange.com/a/87512/50076

Note the code is MIT licensed, which requires attribution.

Code

Available at Github. Algorithm outputs are in test_algo.py.

Answer 2

This is more meat-and-potatoes.

An oft-used algorithm in audio is the pitch detector, which of course is attempting to accurately determine the instantaneous pitch of a note. The note is a waveform that is quasi-periodic.

Referring to the linked answer, the input audio is $x[n]$ and first we need to compute the ASDF (Average Squared Difference Function) which looks a lot like the more commonly referenced AMDF.

$$ Q_x[k, n_0] \triangleq \frac{1}{N} \sum\limits_{n=0}^{N-1} \left(x[n+n_0-\left\lfloor \tfrac{N+k}{2}\right\rfloor] \ - \ x[n+n_0-\left\lfloor \tfrac{N+k}{2}\right\rfloor + k] \right)^2 $$

$\left\lfloor \cdot \right\rfloor$ is the floor() function and, if $k$ is even then $ \left\lfloor \frac{k}{2}\right\rfloor = \left\lfloor \frac{k+1}{2}\right\rfloor = \frac{k}{2} $.

$N$ should be larger than any expected period length. Best that it's twice the longest expected period.

Then, from the ASDF, we can get the autocorrelation function:

$$ R_x[k,n_0] = R_x[0,n_0] - \tfrac12 Q_x[k, n_0] $$ where $$ R_x[0, n_0] \triangleq \frac{1}{N} \sum\limits_{n=0}^{N-1} \Big(x[n+n_0-\left\lfloor \tfrac{N}{2}\right\rfloor]\Big)^2 $$

We know that

$$ Q_x[k, n_0] \ge 0 \qquad \forall k $$ $$ R_x[k, n_0] \le R_x[0, n_0] \qquad \forall k $$

and we know that $Q_x[k, n_0]$ and $R_x[k, n_0]$ are (virtually) even-symmetry functions.

$$ Q_x[k, n_0] = Q_x[-k, n_0] $$ $$ R_x[k, n_0] = R_x[-k, n_0] $$

at least for even $k$ and $N$.

Then there is an important intermediate measure that we often call the normalized autocorrelation.

$$ r_x[k,n_0] \triangleq \frac{R_x[k,n_0]}{R_x[0,n_0]} = 1 - \frac{Q_x[k,n_0]}{2 R_x[0,n_0]} $$

We know that a maximum occurs at $k=0$ and $r_x[0,n_0]=1$. If $x[n]$ is perfectly periodic around $n=n_0$ then there is a $k$ where $r_x[k,n_0] \approx 1$.

If you pass your audio through a DC blocking filter, then we know that the "DC" component of the auto-correlation is also zero and as much area of $r_x[k,n_0]$ is above zero as below zero. So find the value of $k$ that makes $r_x[k,n_0]$ maximum and is not at or around the peak where $k=0$. The YIN algorithm has a dumb way of doing that. Don't bother with YIN.

So there exists a value of $k$, we'll call $k_z$ where $r_x[k,n_0]$ first crosses over from positive values to negative values. So

$$ r_x[k,n_0] \ge 0 \qquad 0 \le k < k_z $$

and

$$ r_x[k_z,n_0] < 0 $$

Now that you know what $k_z$ is, then look for this maximum:

$$ r_x[k_\text{max},n_0] \ge r_x[k,n_0] \qquad \forall k > k_z $$

Now that $r_x[k_\text{max},n_0] \approx 1$ if the waveform is nearly periodic in the neighborhood $n \approx n_0$ $r_x[k_\text{max},n_0]$ can never exceed $1$ and, if you block out all of the DC, then $r_x[k_\text{max},n_0]$ must exceed $0$. We call that value the "pitch confidence" or the "periodicity" measure of the quasi-periodic waveform.

$$ 0 < r_x[k_\text{max},n_0] \le 1 $$

If $r_x[k_\text{max},n_0]$ is big (close to 1), then it's periodic in the neighborhood of $n_0$, when $r_x[k_\text{max},n_0]$ is close to zero, the waveform is crap. Noise. Transient. Something other than periodic.

That's what I would look for.

Answer 3

One domain where continuous (complex) wavelets (CWT) are especially efficient is when you expect the 1D signals to be piecewise regular. When a signal is $C^\alpha$ by pieces, then the $\alpha$ regularity reflects on the scalogram by specific patterns. Moreover, you can expect a decay in wavelet best $M$-term approximation as $M^{-2\alpha}$. So, provided you can estimate correctly a power-law decay, you could recover the original signal regularity.

This may be interesting in your context. First, if the different speech signals that are glued together don't have the same average, the same amplitude, this will show up more easily in the CWT domain, because it enhances diversity. Second, is the segments are normalized, the above may not hold. But if the segments are different regularity, this may show up, albeit in a more subtle way. This also applies to more stochastic phenomena, like fractional Brownian or multifractal processes, and there is literature for their detection.

Of course, these assumptions are almost never known beforehand, and their estimation may be hampered by sampling, discretization and noise.

Between the theoretical context and the real world, I have noticed that complex continuous wavelets can be great tools to detect "a change" at a certain scale in 1D data. My experience in speech is tiny, yet a little bigger in other acoustic/vibration signals. One key is really to identify "on which range of scales change appears", a second key is "can we infer a decision rule on that".

I definitely cannot tell for the 2nd key. For the first one, a couple of picture illustrating that the proper scale may enlighten a tiny change.

1. The loss in amplitude in the noisy signal (top) can be witnessed in a clear way in the sudden decay in the first wavelet spiral (middle), not in the second (bottom)
2. Plots of different time-scale panels for the aforementioned experiment.

Answer 4

I would point on 2 approaches:

1. Parametric Model
   You build a parametric model for the data and then apply simple change detection methods for the paraments. For instance, you may build a Kalman Filter based model and see when the next sample is outside the likelihood range of the prediction. If it is a classic speech you may use parametric models for speech.

2. Non Parametric Models
   You may use simple change detection methods on some transformation of the data. Energy of a frequency bin of the STFT / Wavelet transform. Energy in a window, etc... You may even extract those features and use Isolation Forest of Local Outlier Factor methods. Both are available on SciKit Learn.

Answer 5

I suggest you to elaborate the signal using a wavelet decomposition. If you find the correct wavelet for your task you wll see in the spectogram a band wich will expose you the cut points. This is not sure but you must try many wavelets tecniques (DWT,MODWT,WAVELETS PACKEST ,ecc..) , many waves (HAAR, Daubechies , Orthogonal , Biorthogonal ) and many levels .