diff --git a/src/remucs/tuning.py b/src/remucs/tuning.py new file mode 100644 index 0000000..e7167d9 --- /dev/null +++ b/src/remucs/tuning.py @@ -0,0 +1,92 @@ +from pathlib import Path +from typing import Tuple +from numpy.typing import ArrayLike, NDArray + +import click +import numpy +import soundfile + +from qdft import QDFT +from qdft.fafe import QFAFE +from qdft.scale import Scale + +from remucs.options import RemucsOptions + + +def findpeaks(x: ArrayLike, n: int) -> NDArray: + + x = numpy.atleast_2d(x) + assert len(x.shape) == 2 + assert x.shape[0] > 0 + assert x.shape[1] > 3 + + a = x[..., 0:-3] + y = x[..., 1:-2] + b = x[..., 2:-1] + + i = (y > a) & (y > b) + j = numpy.argpartition(numpy.negative(y * i), n)[..., :n] + + return j + 1 + + +def analyze(src: Path, opts: RemucsOptions) -> Tuple[NDArray, NDArray]: + + if not opts.quiet: + click.echo(f'Analyzing {src.resolve()}') + + x, sr = soundfile.read(src) + x = numpy.atleast_2d(x).mean(axis=-1) + + oldsize = len(x) + newsize = int(numpy.ceil(oldsize / sr) * sr) + x.resize(newsize) + + scale = Scale(440) + bandwidth = (scale.frequency('A0'), scale.frequency('C#8')) + resolution = 12*4 + + qdft = QDFT(samplerate=sr, bandwidth=bandwidth, resolution=resolution) + fafe = QFAFE(qdft) + + batches = numpy.arange(len(x)).reshape((-1, sr)) + estimates = numpy.full(len(x), 440, float) + weights = numpy.zeros(len(x), float) + + numpeaks = 3 + kernel = 1 # int(100e-3 * sr) + + # use qdft.latencies in the next qdft release + roi = [int(numpy.max(qdft.periods[0] - qdft.offsets)), oldsize - 1] + + for batch in batches: + + dfts = qdft.qdft(x[batch]) + magns = numpy.abs(dfts) + freqs = fafe.hz(dfts) + + i = numpy.arange(len(batch)) + j = findpeaks(magns, numpeaks) + + for n, m in zip(i, j): + + if (batch[n] < roi[0]) or (roi[1] < batch[n]): + continue + + estimate = numpy.median(numpy.roll(estimates, batch[n] + kernel)[:kernel]) \ + if kernel > 1 else estimates[batch[n] - 1] + + a = numpy.round(12 * numpy.log2(freqs[n, m] / estimate)) + b = numpy.power(2, a / 12) + c = numpy.power(2, a / 6) + + estimate = numpy.sum(freqs[n, m] * b) / numpy.sum(c) + + if not numpy.isfinite(estimate): + estimates[batch[n]] = estimates[batch[n] - 1] + continue + + estimates[batch[n]] = estimate + weights[batch[n]] = numpy.prod(magns[n, m]) + + return estimates[:oldsize], weights[:oldsize] diff --git a/src/sandbox/tuning.py b/src/sandbox/tuning.py index afbae9a..65e9d28 100644 --- a/src/sandbox/tuning.py +++ b/src/sandbox/tuning.py @@ -1,39 +1,19 @@ # pylint: disable=import-error -# pylint: disable=too-many-locals -# pylint: disable=too-many-statements -# pylint: disable=fixme + +import pathlib import matplotlib.pyplot as plot import numpy as np -import numpy.lib.stride_tricks as tricks import scipy.signal import soundfile -from qdft import QDFT -from qdft.fafe import QFAFE -from qdft.scale import Scale from synth import synth +from remucs import RemucsOptions +from remucs.tuning import analyze np.set_printoptions(suppress=True) -def findpeaks(x, n): - - if n == 1: - return np.argmax(x, axis=-1) - - a = x[..., 0:-3] - y = x[..., 1:-2] - b = x[..., 2:-1] - - i = (y > a) & (y > b) - j = np.argsort(y * i, axis=-1)[..., ::-1] - - assert np.all(np.argmax(y, axis=-1) == j[..., 0]) - - return j[..., :n] + 1 - - def smooth_savgol(x, seconds, samplerate): return scipy.signal.savgol_filter(x, @@ -42,142 +22,48 @@ def smooth_savgol(x, seconds, samplerate): mode='mirror') -def smooth_polyfit(x, degree): - - i = np.arange(len(x)) - p = np.polyfit(i, x, degree) - y = np.poly1d(p)(i) - - print(f'polyfit degree {degree} coeffs {p}') - - return y - - def main(): - cp = 440 - test = f'test.{cp}.wav' - synth(test, a4=cp, tenuto=1) - - samples, samplerate = soundfile.read(test) - - samples = np.mean(samples, axis=-1) \ - if len(np.shape(samples)) > 1 \ - else np.asarray(samples) - - print(f'old length {len(samples)} {len(samples)/samplerate}s') - length = int(np.ceil(samples.size / samplerate) * samplerate) - samples.resize(length) - print(f'new length {len(samples)} {len(samples)/samplerate}s') + cp0 = 440 + test = pathlib.Path(f'test.{cp0}.wav') + synth(test, a4=cp0, tenuto=1) + sr = soundfile.info(test).samplerate - scale = Scale() - bw = (scale.frequency('A0'), scale.frequency('C#8')) - qdft = QDFT(samplerate=samplerate, bandwidth=bw, resolution=12*4) + cp2, weights = analyze(test, RemucsOptions()) - chunks = tricks.sliding_window_view(samples, samplerate)[::samplerate] - data = np.empty((0, qdft.size)) + weights /= np.max(np.abs(weights)) - for i, chunk in enumerate(chunks): - - if not i: - print('0%') - - data = np.vstack((data, qdft.qdft(chunk))) - - print(f'{int(100 * (i + 1) / len(chunks))}%') - - fafe = QFAFE(qdft) - mags = np.abs(data) - freqs = fafe.hz(data) - data = mags + 1j * freqs - - # TODO use qdft.latencies in the next qdft release - latency = int(np.max(qdft.periods[0] - qdft.offsets)) - print(f'max. qdft latency {latency}') - - print(f'old shape {data.shape}') - data = data[latency:-samplerate] - print(f'new shape {data.shape}') - - cp0 = scale.concertpitch - cp1 = np.full(len(data), cp0, float) - - r = np.real(data) - f = np.imag(data) - - i = np.arange(len(data)) - j = findpeaks(r, 3) - k = 1 # int(100e-3 * samplerate) - - for n, m in zip(i, j): - - e = np.median(np.roll(cp1, k)[:k]) \ - if k > 1 else cp1[n-1] - - s = np.round(12 * np.log2(f[n, m] / e)) - - a = f[n, m] - b = np.power(2, s / 12) - c = np.power(2, s / 6) + values = np.round(cp2).astype(int) + minmax = np.min(values), np.max(values) + bins = np.arange(minmax[0], minmax[1] + 1) + edges = np.arange(minmax[0], minmax[1] + 2) - 0.5 + hist = np.histogram(values, bins=edges, weights=weights) + cp1 = bins[np.argmax(hist[0])] + assert hist[0].shape == bins.shape + assert hist[1].shape == edges.shape - cp1[n] = np.sum(a * b) / np.sum(c) - cp1[n] = cp1[n-1] if np.isnan(cp1[n]) else cp1[n] + cp3 = smooth_savgol(cp2, 100e-3, sr) - cp2 = smooth_savgol(cp1, 100e-3, samplerate) - cp3 = smooth_polyfit(cp2, 1) + print(f'cp orig {cp0} est {cp1}') - vals = np.round(cp1).astype(int) - vmin = np.min(vals) - vmax = np.max(vals) - bins = np.arange(vmin, vmax + 1) - edge = np.arange(vmin, vmax + 2) - 0.5 - weig = np.prod(r[i[..., None], j], axis=-1) - hist = np.histogram(vals, bins=edge, weights=weig) - hmax = np.argmax(hist[0]) - bmax = bins[hmax] - assert hist[0].shape == bins.shape - assert hist[1].shape == edge.shape - - res = np.round([ - cp1[0], - cp1[-1], - np.min(cp1), - np.max(cp1), - np.mean(cp1), - np.median(cp1), - bmax - ]).astype(int) - - lab = [ - 'first', - 'last', - 'min', - 'max', - 'avg', - 'med', - 'hist', - ] - - print('\n'.join(map(str, zip(lab, res)))) - - plot.figure(test + ' hist') - plot.hist(vals, bins=edge) - plot.axvline(bmax, label='bmax', linestyle='--') - plot.axvline(cp, label='cp', linestyle='--', color='m') + plot.figure(f'{test} hist') + plot.hist(values, bins=edges, weights=weights) + plot.axvline(cp0, label='cp0', linestyle='-', color='c') + plot.axvline(cp1, label='cp1', linestyle='--', color='m') plot.legend() plot.tight_layout() - plot.figure(test) - plot.plot(cp1, label='cp1') - plot.plot(cp2, label='cp2', linestyle='--') - plot.plot(cp3, label='cp3', linestyle='--') - plot.axhline(cp, label='cp', linestyle='--', color='m') + plot.figure(f'{test} plot') + plot.plot(cp2) + plot.plot(cp3, linestyle='--') + plot.axhline(cp0, label='cp0', linestyle='-', color='c') + plot.axhline(cp1, label='cp1', linestyle='--', color='m') plot.legend() plot.tight_layout() plot.show() - assert res[-1] == cp + assert cp1 == cp0 if __name__ == '__main__':