Merge 9fdad73 into 397c724

librosa/core/__init__.py

@@ -30,6 +30,8 @@
     pseudo_cqt
     fmt
 
+    harmonics
+
     phase_vocoder
 
     magphase

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.
 

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)])