loop.py

import sys
import numpy
from tqdm import tqdm
from truncate import read_wave_file
from audio_helpers import fundamental_frequency

QUANTIZE_FACTOR = 8


def compare_windows(window_a, window_b):
    return numpy.sqrt(numpy.mean(numpy.power(window_a - window_b, 2)))


def slide_window(file, period, start_at=0, end_before=0):
    for power in reversed(xrange(7, 10)):
        multiple = 2 ** power
        window_size = int(period * multiple)
        # Uncomment this to search from the start_at value to the end_before
        # rather than just through one window's length
        # end_range = len(file) - (window_size * 2) - end_before
        end_range = start_at + window_size
        for i in xrange(start_at, end_range):
            yield power, i, window_size


def window_match(file):
    period = (1.0 / fundamental_frequency(file, 1)) * 2
    print period, 'period in samples'

    winner = None

    window_positions = list(
        slide_window(file, period, len(file) / 2, len(file) / 8)
    )
    for power, i, window_size in tqdm(window_positions):
        window_start = find_similar_sample_index(file, i, i + window_size)
        window_end = find_similar_sample_index(file, i, i + (window_size * 2))
        effective_size = window_end - window_start

        difference = compare_windows(
            file[i:i + effective_size],
            file[window_start:window_end]
        ) / effective_size
        if winner is None or difference < winner[0]:
            winner = (
                difference,
                effective_size,
                i,
                abs(file[i] - file[window_start])
            )
            print 'new winner', winner

    lowest_difference, winning_window_size, winning_index, gap = winner

    print "Best loop match:", lowest_difference
    print "window size", winning_window_size
    print "winning index", winning_index
    print "winning gap", gap
    return winning_index, winning_window_size


def slope_at_index(file, i):
    return (file[i + 1] - file[i - 1]) / 2


def find_similar_sample_index(
    file,
    reference_index,
    search_around_index,
    search_size=100  # samples
):
    reference_slope = slope_at_index(file, reference_index) > 0
    best_match = None
    search_range = xrange(
        search_around_index - search_size,
        search_around_index + search_size
    )
    for i in search_range:
        slope = slope_at_index(file, i) > 0
        if slope != reference_slope:
            continue

        abs_diff = abs(file[i] - file[reference_index])

        if best_match is not None:
            _, best_abs_diff = best_match
            if abs_diff < best_abs_diff:
                best_match = (i, abs_diff)
        else:
            best_match = (i, abs_diff)
    return best_match[0] if best_match is not None else search_around_index


def zero_crossing_match(file):
    period = (1.0 / fundamental_frequency(file, 1)) * 2
    print period, 'period in samples'

    period_multiple = 64
    period = period * period_multiple

    for i in reversed(xrange(2 * len(file) / 3, 5 * len(file) / 6)):
        if file[i] >= 0 and file[i + 1] < 0 and \
                file[int(i + period)] >= 0 and \
                file[int(i + 1 + period)] < 0 and \
                file[int(i + period * 2)] >= 0 and \
                file[int(i + 1 + period * 2)] < 0:
            return i, int(period)


def fast_autocorrelate(x):
    """
    Compute the autocorrelation of the signal, based on the properties of the
    power spectral density of the signal.

    Note that the input's length may be reduced before the correlation is
    performed due to a pathalogical case in numpy.fft:
    http://stackoverflow.com/a/23531074/679081

    > The FFT algorithm used in np.fft performs very well (meaning O(n log n))
    > when the input length has many small prime factors, and very bad
    > (meaning a naive DFT requiring O(n^2)) when the input size is a prime
    > number.
    """

    # This is one simple way to ensure that the input array
    # has a length with many small prime factors, although it
    # doesn't guarantee that (also hopefully we don't chop too much)
    optimal_input_length = int(numpy.sqrt(len(x))) ** 2
    x = x[:optimal_input_length]
    xp = x - numpy.mean(x)
    f = numpy.fft.fft(xp)
    p = numpy.absolute(numpy.power(f, 2))
    pi = numpy.fft.ifft(p)
    result = numpy.real(pi)[:x.size / 2] / numpy.sum(numpy.power(xp, 2))
    return result


def find_argmax_after(file, offset):
    return numpy.argmax(file[offset:]) + offset


def autocorrelated_loop(file, search_start, min_loop_width_in_seconds=0.2):
    # Strategy:
    #  1) run an autocorrelation on the file.
    #  3) Find argmax of the autocorrelation
    #  4) define some peak_width and find the next highest peak after current
    #  5) define the loop bounds as from the first peak to the second peak
    #  6) massage the loop bounds using find_similar_sample_index
    #  7) ???
    #  8) Profit!
    autocorrelation = fast_autocorrelate(file)
    return find_loop_from_autocorrelation(
        file,
        autocorrelation,
        search_start,
        min_loop_width_in_seconds
    )


def find_loop_from_autocorrelation(
    file,
    autocorrelation,
    search_start,
    min_loop_width_in_seconds=0.2,
    sample_rate=48000
):
    search_start /= 2
    max_autocorrelation_peak_width = int(
        min_loop_width_in_seconds * sample_rate
    )
    loop_start = find_argmax_after(autocorrelation, search_start)
    loop_end = find_argmax_after(
        autocorrelation,
        loop_start + max_autocorrelation_peak_width
    )

    loop_end = find_similar_sample_index(file, loop_start, loop_end) - 1
    return loop_start, (loop_end - loop_start)


def minimize(iterable, callable):
    best_result = None
    best_score = None
    for x in iterable:
        if x:
            score = callable(*x)
            if best_score is None or score < best_score:
                best_score = score
                best_result = x
    return best_result


def autocorrelate_loops(file, sample_rate):
    autocorrelation = fast_autocorrelate(file)
    search_points = [
        3 * len(file) / 4,
        2 * len(file) / 3,
        len(file) / 2,
        len(file) / 3,
    ]
    loop_widths = [0.2, 0.4, 0.6, 0.8, 1.0, 1.5, 2, 2.5, 3.]
    for search_point in search_points:
        for width in loop_widths:
            try:
                yield find_loop_from_autocorrelation(
                    file, autocorrelation,
                    search_point, width, sample_rate)
            except ValueError:
                # We couldn't search for a loop width of that size.
                pass
    yield None


def find_loop_points(data, sample_rate):
    channel = data[0]

    result = minimize(
        autocorrelate_loops(channel, sample_rate),
        lambda start, length: abs(channel[start] - channel[start + length])
    )

    if result:
        loop_start, loop_size = result
        return loop_start, loop_start + loop_size


def process(aif, sample_rate=48000):
    file = read_wave_file(aif)

    # loop_start, loop_size = window_match(file)
    # loop_start, loop_size = zero_crossing_match(file)
    loop_start, loop_end = find_loop_points(file)
    loop_size = loop_end - loop_start

    file = file[0]

    print 'start, end', loop_start, loop_end

    plt.plot(file[loop_start:loop_end])
    plt.plot(file[loop_end:loop_start + (2 * loop_size)])
    plt.show()

    plt.plot(file[
        loop_start - (sample_rate * 2):
        loop_start + (sample_rate * 2)
    ])
    plt.axvline(sample_rate * 2)
    plt.axvline((sample_rate * 2) + loop_size)
    plt.show()

if __name__ == "__main__":
    import matplotlib.pyplot as plt
    process(sys.argv[1])