Ripple sounds

ripple.py

# -*- coding: utf-8 -*-

from warnings import warn
import numpy as np
from scipy.signal import butter, lfilter
from scipy.stats import norm, uniform
import matplotlib.pyplot as plt
import brian
import brian.hears as hears


def butter_lowpass(cutoff, fs=44100., order=4):
    """Designs a low-pass Butterworth filter."""
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a


def butter_lowpass_filter(data, cutoff, fs, order=4):
    """Applies a low-pass Butterworth filter to the data."""
    b, a = butter_lowpass(cutoff, fs, order)
    y = lfilter(b, a, data)
    return y


def random_walk(duration, cutoff, bounds, fs=44100., order=4):
    """Calculates the cumulative sum of duration * fs random gaussian draws,
    low-pass filters the vector and scales this between bounds. The result is
    a smooth random walk."""
    t = np.linspace(0, 1 * duration, duration * fs)
    y = norm.rvs(size=t.size)
    y = butter_lowpass_filter(y, cutoff, fs, order)
    w = np.cumsum(y)
    w = w - w.min()
    w = w / w.max()
    w = w * (bounds[1] - bounds[0]) + bounds[0]
    return w


def ripple_sound(duration, ntones, omega, w, randomise_gamma=True, d=1.,
                 ripple_phi=0., fmin=250., fmax=8000., fs=44100.,
                 filename=None, return_env=False):
    """Synthesises a ripple sound. Requires the duration of the sound in s,
    the number of tones to equally space between fmin and fmax, omega and w
    values. Omega and w can be scalars or vectors on length duration * fs.
    Returns the synthesised sound, and optionally the envelope vector. If a
    filename is given, the sound is saved at 90% max amplitude."""

    if ntones > 1000:
        warn('ntones exceeds 1000; this may take a very long time')

    if randomise_gamma:
        gamma = uniform.rvs(0, 1, size=ntones)
    else:
        gamma = np.ones(ntones)

    if hasattr(omega, '__iter__'):
        if len(omega) != duration * fs:
            raise Exception('omega vector has incorrect length')
    else:
        omega = np.tile(omega, duration * fs)

    if hasattr(w, '__iter__'):
        if len(w) != duration * fs:
            raise Exception('w vector has incorrect length')
    else:
        w = np.tile(w, duration * fs)

    t = np.linspace(0, duration, duration * fs)
    i = np.arange(ntones) + 1
    f = fmin * (fmax / fmin)**((i - 1)/(float(ntones) - 1))
    phi = uniform.rvs(0, 2 * np.pi, size=ntones)

    T = np.tile(t, (ntones, 1))
    F = np.tile(f, (t.size, 1)).T
    Gamma = np.tile(gamma, (t.size, 1)).T
    Phi = np.tile(phi, (t.size, 1)).T
    Omega = np.tile(omega, (ntones, 1))
    wprime = np.cumsum(w) / float(fs)
    WP = np.tile(wprime, (ntones, 1))
    X = np.log2(F / fmin)

    A = 1 + d * np.sin(2 * np.pi * (WP * T + Omega * X) + ripple_phi)
    S = Gamma * np.sin(2 * np.pi * F * T + Phi) / np.sqrt(F)
    Y = A * S
    y = np.sum(Y, axis=0)
    y = y / np.abs(y).max()

    if filename:
        y = y * .9
        sound = hears.Sound(y)
        sound.save(filename)

    if return_env:
        return y, A
    else:
        return y


def main():
    """Creates a single dynamic moving ripple sound, shows its envelope and
    spectrogram, and saves."""
    duration = 2.
    ntones = 1000
    omega = random_walk(duration, 2, (0, 4.5))
    w = random_walk(duration, 2, (-6, 6))
    y, a = ripple_sound(duration, ntones, omega, w, filename='ripple_4.wav',
                        return_env=True)
    plt.pcolormesh(a)
    plt.xlabel('Time (s)')
    plt.ylabel('$i$')
    plt.xlim(0, duration * 44100.)
    ticks = np.array([0.0, 0.4, 0.8, 1.2, 1.6, 1.8])
    plt.xticks(ticks*44100, ticks)
    plt.show()
    y = y * .9
    sound = hears.Sound(y)
    sound.spectrogram(0, 10000)
    plt.show()


if __name__ == '__main__':
    main()