Merge eee3975 into fe7f444

librosa/effects.py

@@ -27,18 +27,21 @@
     :toctree: generated/
 
     remix
+    trim
+    split
 """
 
 import numpy as np
 
 from . import core
 from . import decompose
+from . import feature
 from . import util
 from .util.exceptions import ParameterError
 
 __all__ = ['hpss', 'harmonic', 'percussive',
            'time_stretch', 'pitch_shift',
-           'remix']
+           'remix', 'trim', 'split']
 
 
 def hpss(y, **kwargs):
@@ -77,7 +80,7 @@ def hpss(y, **kwargs):
     >>> y, sr = librosa.load(librosa.util.example_audio_file())
     >>> y_harmonic, y_percussive = librosa.effects.hpss(y)
 
-    >>> # Get a more isolated percussive component by widening its margin 
+    >>> # Get a more isolated percussive component by widening its margin
     >>> y_harmonic, y_percussive = librosa.effects.hpss(y, margin=(1.0,5.0))
 
     '''
@@ -94,6 +97,7 @@ def hpss(y, **kwargs):
 
     return y_harm, y_perc
 
+
 def harmonic(y, **kwargs):
     '''Extract harmonic elements from an audio time-series.
 
@@ -117,14 +121,13 @@ def harmonic(y, **kwargs):
 
     Examples
     --------
-    >>> # Extract harmonic component 
+    >>> # Extract harmonic component
     >>> y, sr = librosa.load(librosa.util.example_audio_file())
     >>> y_harmonic = librosa.effects.harmonic(y)
 
     >>> # Use a margin > 1.0 for greater harmonic separation
     >>> y_harmonic = librosa.effects.harmonic(y, margin=3.0)
 
-
     '''
 
     # Compute the STFT matrix
@@ -138,6 +141,7 @@ def harmonic(y, **kwargs):
 
     return y_harm
 
+
 def percussive(y, **kwargs):
     '''Extract percussive elements from an audio time-series.
 
@@ -160,15 +164,14 @@ def percussive(y, **kwargs):
     librosa.decompose.hpss : HPSS for spectrograms
 
     Examples
-    --------    
-    >>> # Extract percussive component 
+    --------
+    >>> # Extract percussive component
     >>> y, sr = librosa.load(librosa.util.example_audio_file())
     >>> y_percussive = librosa.effects.percussive(y)
 
     >>> # Use a margin > 1.0 for greater percussive separation
     >>> y_percussive = librosa.effects.percussive(y, margin=3.0)
 
-
     '''
 
     # Compute the STFT matrix
@@ -374,3 +377,165 @@ def remix(y, intervals, align_zeros=True):
         y_out.append(y[clip])
 
     return np.concatenate(y_out, axis=-1)
+
+
+def _signal_to_frame_nonsilent(y, frame_length=2048, hop_length=512, top_db=60,
+                               ref_power=np.max):
+    '''Frame-wise non-silent indicator for audio input.
+
+    This is a helper function for `trim` and `split`.
+
+    Parameters
+    ----------
+    y : np.ndarray, shape=(n,) or (2,n)
+        Audio signal, mono or stereo
+
+    frame_length : int > 0
+        The number of samples per frame
+
+    hop_length : int > 0
+        The number of samples between frames
+
+    top_db : number > 0
+        The threshold (in decibels) below reference to consider as
+        silence
+
+    ref_power : callable or float
+        The reference power
+
+    Returns
+    -------
+    non_silent : np.ndarray, shape=(m,), dtype=bool
+        Indicator of non-silent frames
+    '''
+    # Convert to mono
+    y_mono = core.to_mono(y)
+
+    # Compute the MSE for the signal
+    mse = feature.rmse(y=y_mono, n_fft=frame_length, hop_length=hop_length)**2
+
+    return (core.logamplitude(mse.squeeze(),
+                              ref_power=ref_power,
+                              top_db=None) > - top_db)
+
+
+def trim(y, top_db=60, ref_power=np.max, frame_length=2048, hop_length=512):
+    '''Trim leading and trailing silence from an audio signal.
+
+    Parameters
+    ----------
+    y : np.ndarray, shape=(n,) or (2,n)
+        Audio signal, can be mono or stereo
+
+    top_db : number > 0
+        The threshold (in decibels) below reference to consider as
+        silence
+
+    ref_power : number or callable
+        The reference power.  By default, it uses `np.max` and compares
+        to the peak power in the signal.
+
+    frame_length : int > 0
+        The number of samples per analysis frame
+
+    hop_length : int > 0
+        The number of samples between analysis frames
+
+    Returns
+    -------
+    y_trimmed : np.ndarray, shape=(m,) or (2, m)
+        The trimmed signal
+
+    index : np.ndarray, shape=(2,)
+        the slice of `y` corresponding to the non-silent region:
+        `y_trimmed = y[index]`.
+
+    Examples
+    --------
+    >>> # Load some audio
+    >>> y, sr = librosa.load(librosa.util.example_audio_file())
+    >>> # Trim the beginning and ending silence
+    >>> yt, sl = librosa.effects.trim(y)
+    >>> # Print the durations
+    >>> print(librosa.get_duration(y), librosa.get_duration(yt))
+    61.45886621315193 60.58086167800454
+    '''
+
+    non_silent = _signal_to_frame_nonsilent(y,
+                                            frame_length=frame_length,
+                                            hop_length=hop_length,
+                                            ref_power=ref_power,
+                                            top_db=top_db)
+
+    nonzero = np.flatnonzero(non_silent)
+
+    # Compute the start and end positions
+    # End position goes one frame past the last non-zero
+    start = int(core.frames_to_samples(nonzero[0], hop_length))
+    end = min(y.shape[-1],
+              int(core.frames_to_samples(nonzero[-1] + 1, hop_length)))
+
+    # Build the mono/stereo index
+    full_index = [slice(None)] * y.ndim
+    full_index[-1] = slice(start, end)
+
+    return y[full_index], np.asarray([start, end])
+
+
+def split(y, top_db=60, ref_power=np.max, frame_length=2048, hop_length=512):
+    '''Split an audio signal into non-silent intervals.
+
+    Parameters
+    ----------
+    y : np.ndarray, shape=(n,) or (2, n)
+        An audio signal
+
+    top_db : number > 0
+        The threshold (in decibels) below reference to consider as
+        silence
+
+    ref_power : number or callable
+        The reference power.  By default, it uses `np.max` and compares
+        to the peak power in the signal.
+
+    frame_length : int > 0
+        The number of samples per analysis frame
+
+    hop_length : int > 0
+        The number of samples between analysis frames
+
+    Returns
+    -------
+    intervals : np.ndarray, shape=(m, 2)
+        `intervals[i] == (start_i, end_i)` are the start and end time
+        (in samples) if the `i`th non-silent interval.
+    '''
+
+    non_silent = _signal_to_frame_nonsilent(y,
+                                            frame_length=frame_length,
+                                            hop_length=hop_length,
+                                            ref_power=ref_power,
+                                            top_db=top_db)
+
+    # Interval slicing, adapted from
+    # https://stackoverflow.com/questions/2619413/efficiently-finding-the-interval-with-non-zeros-in-scipy-numpy-in-python
+    # Find points where the sign flips
+    edges = np.flatnonzero(np.diff(non_silent.astype(int)))
+
+    # Pad back the sample lost in the diff
+    edges = [edges + 1]
+
+    # If the first frame had high energy, count it
+    if non_silent[0]:
+        edges.insert(0, [0])
+
+    # Likewise for the last frame
+    if non_silent[-1]:
+        edges.append([len(non_silent)])
+
+    # Convert from frames to samples
+    edges = core.frames_to_samples(np.concatenate(edges),
+                                   hop_length=hop_length)
+
+    # Stack the results back as an ndarray
+    return edges.reshape((-1, 2))

tests/test_effects.py

@@ -139,3 +139,86 @@ def test_harmonic():
     yh2 = librosa.effects.harmonic(y)
 
     assert np.allclose(yh1, yh2)
+
+
+def test_trim():
+
+    def __test(y, top_db, ref_power, trim_duration):
+        yt, idx = librosa.effects.trim(y, top_db=top_db,
+                                       ref_power=ref_power)
+
+        # Test for index position
+        fidx = [slice(None)] * y.ndim
+        fidx[-1] = slice(*idx.tolist())
+        assert np.allclose(yt, y[fidx])
+
+        # Verify logamp
+        rms = librosa.feature.rmse(librosa.to_mono(yt))
+        logamp = librosa.logamplitude(rms**2, ref_power=ref_power, top_db=None)
+        assert np.all(logamp > - top_db)
+
+        # Verify logamp
+        rms_all = librosa.feature.rmse(librosa.to_mono(y)).squeeze()
+        logamp_all = librosa.logamplitude(rms_all**2, ref_power=ref_power,
+                                          top_db=None)
+
+        start = int(librosa.samples_to_frames(idx[0]))
+        stop = int(librosa.samples_to_frames(idx[1]))
+        assert np.all(logamp_all[:start] <= - top_db)
+        assert np.all(logamp_all[stop:] <= - top_db)
+
+        # Verify duration
+        duration = librosa.get_duration(yt)
+        assert np.allclose(duration, trim_duration, atol=1e-1), duration
+
+    # construct 5 seconds of stereo silence
+    # Stick a sine wave in the middle three seconds
+    sr = float(22050)
+    trim_duration = 3.0
+    y = np.sin(2 * np.pi * 440. * np.arange(0, trim_duration * sr) / sr)
+    y = librosa.util.pad_center(y, 5 * sr)
+    y = np.vstack([y, np.zeros_like(y)])
+
+    for top_db in [60, 40, 20]:
+        for ref_power in [1, np.max]:
+            # Test stereo
+            yield __test, y, top_db, ref_power, trim_duration
+            # Test mono
+            yield __test, y[0], top_db, ref_power, trim_duration
+
+
+def test_split():
+
+    def __test(hop_length, frame_length, top_db):
+
+        intervals = librosa.effects.split(y,
+                                          top_db=top_db,
+                                          frame_length=frame_length,
+                                          hop_length=hop_length)
+
+        int_match = librosa.util.match_intervals(intervals, idx_true)
+
+        for i in range(len(intervals)):
+            i_true = idx_true[int_match[i]]
+
+            assert np.all(np.abs(i_true - intervals[i]) <= frame_length), intervals[i]
+
+    # Make some high-frequency noise
+    sr = 8192
+
+    y = np.ones(10 * sr)
+    y[::2] *= -1
+
+    # Zero out all but two intervals
+    y[:sr] = 0
+    y[2 * sr:3 * sr] = 0
+    y[4 * sr:] = 0
+
+    # The true non-silent intervals
+    idx_true = np.asarray([[sr, 2 * sr],
+                           [3 * sr, 4 * sr]])
+
+    for frame_length in [1024, 2048, 4096]:
+        for hop_length in [256, 512, 1024]:
+            for top_db in [20, 60, 80]:
+                yield __test, hop_length, frame_length, top_db