Merge 05dd0c7 into fb87932

librosa/core/spectrum.py

@@ -9,14 +9,16 @@
 import scipy.interpolate
 import six
 
+from . import harmonic
 from . import time_frequency
+from .. import decompose
 from .. import cache
 from .. import util
 from ..util.exceptions import ParameterError
 from ..filters import get_window
 
 __all__ = ['stft', 'istft', 'magphase',
-           'ifgram',
+           'ifgram', 'salience',
            'phase_vocoder',
            'logamplitude', 'perceptual_weighting',
            'fmt']
@@ -54,7 +56,7 @@ def stft(y, n_fft=2048, hop_length=None, win_length=None, window='hann',
         If unspecified, defaults to ``win_length = n_fft``.
 
     window : string, tuple, number, function, or np.ndarray [shape=(n_fft,)]
-        - a window specification (string, tuple, or number); 
+        - a window specification (string, tuple, or number);
           see `scipy.signal.get_window`
         - a window function, such as `scipy.signal.hanning`
         - a vector or array of length `n_fft`
@@ -439,6 +441,90 @@ def ifgram(y, sr=22050, n_fft=2048, hop_length=None, win_length=None,
     return if_gram, stft_matrix
 
 
+def salience(S, freqs, harmonics, weights, peak_filter=True, use_hpss=True):
+    """Harmonic salience function.
+
+    Parameters
+    ----------
+    S : np.ndarray [shape=(d, n)]
+        input time frequency representation (stft, ifgram, etc).
+        May be real (magnitude) or complex.
+
+    freqs : np.ndarray, shape=(X.shape[axis])
+        The frequency values corresponding to X's elements along the
+        chosen axis.
+
+    harmonics : list-like, non-negative
+        Harmonics to include in salience computation.  The first harmonic (1)
+        corresponds to `S` itself. Values less than one (e.g., 1/2) correspond
+        to sub-harmonics.
+
+    weights : list-like
+        The weight to apply to each harmonic in the summation.
+        Must be the same length as `harmonics`.
+
+    peak_filter : bool
+        If true, computes harmonic summation only on peaks in frequency.
+        Otherwise computes harmonic summation over the full spectrum.
+        Defaults to True.
+
+    use_hpss : bool
+        If true, the salience function is computed on the harmonic portion
+        of the input representation. Defaults to True.
+
+    Returns
+    -------
+    S_sal : np.ndarray, shape=(len(h_range), [x.shape])
+        `S_sal` will have the same shape as `S`, and measure
+        the overal harmonic energy at each frequency.
+
+
+    See Also
+    --------
+    core.harmonics
+
+
+    Examples
+    --------
+    >>> y, sr = librosa.load(librosa.util.example_audio_file(),
+    ...                      duration=15, offset=30)
+    >>> S = np.abs(librosa.stft(y))
+    >>> freqs = librosa.core.fft_frequencies(sr)
+    >>> harmonics = [1./3, 1./2, 1, 2, 3, 4]
+    >>> weights = [-0.5, -1.0, 1.0, 0.5, 0.33, 0.25]
+    >>> S_sal = librosa.salience(S, freqs, harmonics, weights)
+    >>> print(S_sal.shape)
+    (1025, 646)
+
+    >>> plt.figure()
+    >>> librosa.display.specshow(librosa.logamplitude(S_sal**2,
+    ...                                               ref_power=S_sal.max()*2),
+    ...                          sr=sr, y_axis='log')
+    >>> plt.tight_layout()
+
+    """
+
+    if np.iscomplexobj(S):
+        S = np.abs(S)
+
+    if use_hpss:
+        S, _ = decompose.hpss(S)
+
+    weights = np.array(weights, dtype=float)
+
+    S_harm = harmonic.harmonics(np.abs(S), freqs, harmonics)
+    S_peaks = scipy.signal.argrelmax(np.abs(S), axis=0)
+
+    if peak_filter:
+        peak_mask = np.ones(S_harm.shape)
+        peak_mask[:, S_peaks[0], S_peaks[1]] = 0
+        S_harm[peak_mask.astype(bool)] = 0
+
+    S_sal = (S_harm.T * weights).T.sum(axis=0)
+
+    return S_sal
+
+
 def magphase(D):
     """Separate a complex-valued spectrogram D into its magnitude (S)
     and phase (P) components, so that `D = S * P`.

tests/test_core.py

@@ -279,6 +279,16 @@ def __test(ref_power, clip):
             yield tf, ref_power, clip
 
 
+def test_salience():
+    (y, sr) = librosa.load('data/test1_22050.wav')
+    S = librosa.stft(y)
+    freqs = librosa.core.fft_frequencies(sr)
+    harmonics = [1, 2]
+    weights = [1.0, 0.5]
+    S_sal = librosa.core.salience(S, freqs, harmonics, weights)
+    #TODO
+
+
 def test_magphase():
 
     (y, sr) = librosa.load('data/test1_22050.wav')