From 2eae412d4ce5eff6bebaa82de28e09438ebbcf77 Mon Sep 17 00:00:00 2001 From: Brian McFee Date: Wed, 17 Aug 2016 15:35:17 -0400 Subject: [PATCH 1/2] added harmonic calculator. needs tests --- librosa/core/__init__.py | 2 + librosa/core/spectrum.py | 119 ++++++++++++++++++++++++++++++++++++++- 2 files changed, 120 insertions(+), 1 deletion(-) diff --git a/librosa/core/__init__.py b/librosa/core/__init__.py index 5e74d7001e..6e0907a302 100644 --- a/librosa/core/__init__.py +++ b/librosa/core/__init__.py @@ -30,6 +30,8 @@ pseudo_cqt fmt + harmonics + phase_vocoder magphase diff --git a/librosa/core/spectrum.py b/librosa/core/spectrum.py index d00c8d64c3..564905fec7 100644 --- a/librosa/core/spectrum.py +++ b/librosa/core/spectrum.py @@ -18,7 +18,7 @@ 'ifgram', 'phase_vocoder', 'logamplitude', 'perceptual_weighting', - 'fmt'] + 'fmt', 'harmonics'] @cache(level=20) @@ -993,6 +993,123 @@ def fmt(y, t_min=0.5, n_fmt=None, kind='cubic', beta=0.5, over_sample=1, axis=-1 return result[idx] * np.sqrt(n) / n_fmt +def harmonics(X, freqs, h_range, kind='slinear', fill_value=0, axis=0): + '''Built a harmonic tensor from a time-frequency representation. + + Parameters + ---------- + X : np.ndarray, real-valued + The input energy + + freqs : np.ndarray, shape=(X.shape[axis]) + The frequency values corresponding to X's elements along the + chosen axis. + + h_range : list-like, non-negative + Harmonics to compute. The first harmonic (1) corresponds to `X` + itself. + Values less than one (e.g., 1/2) correspond to sub-harmonics. + + kind : str + Interpolation type. See `scipy.interpolate.interp1d`. + + axis : int + The axis along which to compute harmonics + + fill_value : float + The value to fill when extrapolating beyond the observed + frequency range. + + Returns + ------- + X_harm : np.ndarray, shape=(len(h_range), [X.shape]) + `X_harm[i]` will have the same shape as `X`, and measure + the energy at the `h_range[i]` harmonic of each frequency. + + See Also + -------- + scipy.interpolate.interp1d + + + Examples + -------- + Estimate the harmonics of a time-averaged tempogram + + >>> y, sr = librosa.load(librosa.util.example_audio_file(), + ... duration=15, offset=30) + >>> # Compute the time-varying tempogram and average over time + >>> tempi = np.mean(librosa.feature.tempogram(y=y, sr=sr), axis=1) + >>> # We'll measure the first five harmonics + >>> h_range = [1, 2, 3, 4, 5] + >>> f_tempo = librosa.tempo_frequencies(len(tempi), sr=sr) + >>> # Build the harmonic tensor + >>> t_harmonics = librosa.harmonics(tempi, f_tempo, h_range) + >>> # We'll ignore any nans or infs + >>> t_harmonics[np.isnan(t_harmonics)] = 0 + >>> print(t_harmonics.shape) + (5, 384) + + >>> # And plot the results + >>> import matplotlib.pyplot as plt + >>> plt.figure() + >>> librosa.display.specshow(t_harmonics, x_axis='tempo', sr=sr) + >>> plt.yticks(0.5 + np.arange(len(h_range)), + ... ['{:.3g}'.format(_) for _ in h_range]) + >>> plt.ylabel('Harmonic') + >>> plt.xlabel('Tempo (BPM)') + >>> plt.tight_layout() + + We can also compute frequency harmonics for spectrograms. + To calculate subharmonic energy, use values < 1. + + >>> h_range = [1./3, 1./2, 1, 2, 3, 4] + >>> S = np.abs(librosa.stft(y)) + >>> fft_freqs = librosa.fft_frequencies(sr=sr) + >>> S_harm = librosa.harmonics(S, fft_freqs, h_range, axis=0) + >>> print(S_harm.shape) + (6, 1025, 646) + + >>> plt.figure() + >>> for i, _sh in enumerate(S_harm, 1): + ... plt.subplot(3,2,i) + ... librosa.display.specshow(librosa.logamplitude(_sh**2, + ... ref_power=S.max()**2), + ... sr=sr, y_axis='log') + ... plt.title('h={:.3g}'.format(h_range[i-1])) + ... plt.yticks([]) + >>> plt.tight_layout() + ''' + + # Note: this only works for fixed-grid, 1d interpolation + f_interp = scipy.interpolate.interp1d(freqs, X, + kind=kind, + axis=axis, + bounds_error=False, + fill_value=fill_value) + + # X_out will be the same shape as X, plus a leading + # axis that has length = len(h_range) + out_shape = [len(h_range)] + out_shape.extend(X.shape) + harmonic_out = np.zeros(out_shape, dtype=X.dtype) + + idx_out = [slice(None)] * harmonic_out.ndim + + # Iterate over the harmonics range + for h_index, harmonic in enumerate(h_range): + idx_out[0] = h_index + + # Iterate over frequencies + for f_index, frequency in enumerate(freqs): + # Offset the output axis by 1 to account for the harmonic index + idx_out[1 + axis] = f_index + + # Estimate the harmonic energy at this frequency across time + harmonic_out[tuple(idx_out)] = f_interp(harmonic * frequency) + + return harmonic_out + + def _spectrogram(y=None, S=None, n_fft=2048, hop_length=512, power=1): '''Helper function to retrieve a magnitude spectrogram. From 9fdad735ecabc46659bda07bfc1502a263c25aa8 Mon Sep 17 00:00:00 2001 From: Brian McFee Date: Wed, 17 Aug 2016 15:58:34 -0400 Subject: [PATCH 2/2] added tests for harmonics --- tests/test_core.py | 48 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/tests/test_core.py b/tests/test_core.py index 1db9097b1c..f0a2ee6b8e 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -852,3 +852,51 @@ def test_fmt_axis(): f2 = librosa.fmt(y.T, axis=0).T assert np.allclose(f1, f2) + + +def test_harmonics_1d(): + + x = np.arange(16) + y = np.linspace(-8, 8, num=len(x), endpoint=False)**2 + + h = [0.25, 0.5, 1, 2, 4] + + yh = librosa.harmonics(y, x, h) + + eq_(yh.shape[1:], y.shape) + eq_(yh.shape[0], len(h)) + for i in range(len(h)): + if h[i] <= 1: + # Check that subharmonics match + step = int(1./h[i]) + vals = yh[i, ::step] + assert np.allclose(vals, y[:len(vals)]) + else: + # Else check that harmonics match + step = h[i] + vals = y[::step] + assert np.allclose(vals, yh[i, :len(vals)]) + + +def test_harmonics_2d(): + + x = np.arange(16) + y = np.linspace(-8, 8, num=len(x), endpoint=False)**2 + y = np.tile(y, (5, 1)).T + h = [0.25, 0.5, 1, 2, 4] + + yh = librosa.harmonics(y, x, h, axis=0) + + eq_(yh.shape[1:], y.shape) + eq_(yh.shape[0], len(h)) + for i in range(len(h)): + if h[i] <= 1: + # Check that subharmonics match + step = int(1./h[i]) + vals = yh[i, ::step] + assert np.allclose(vals, y[:len(vals)]) + else: + # Else check that harmonics match + step = h[i] + vals = y[::step] + assert np.allclose(vals, yh[i, :len(vals)])