Skip to content

Commit

Permalink
Refactor tuning prototype #1
Browse files Browse the repository at this point in the history
  • Loading branch information
jurihock committed Mar 31, 2024
1 parent 57831d9 commit c716d3d
Show file tree
Hide file tree
Showing 2 changed files with 122 additions and 144 deletions.
92 changes: 92 additions & 0 deletions src/remucs/tuning.py
Original file line number Diff line number Diff line change
@@ -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]
174 changes: 30 additions & 144 deletions src/sandbox/tuning.py
Original file line number Diff line number Diff line change
@@ -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,
Expand All @@ -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__':
Expand Down

0 comments on commit c716d3d

Please sign in to comment.