Pitch class and height

import os
import urllib.request

import matplotlib.pyplot as plt
import numpy as np

import sonore as so

plt.rcParams.update({"font.size": 9, "axes.titlesize": 10, "figure.dpi": 100})


def fetch(path):
    """A file from the sonore repository, by its path there: the local copy when this runs from
    the repository root, otherwise downloaded from GitHub to the same relative path."""
    if not os.path.exists(path):
        os.makedirs(os.path.dirname(path) or ".", exist_ok=True)
        urllib.request.urlretrieve("https://raw.githubusercontent.com/choyun1/sonore/main/" + path, path)
    return path


def finish(snd):
    """How every sound in the gallery is played: 5 ms ramps, RMS 0.1, peak at most 0.95."""
    snd = snd.ramp(5e-3).normalize(rms=0.1)
    return snd.normalize(peak=0.95) if snd.peak > 0.95 else snd

FS = 22050
LEAD = 0.5  # silence before every sound [s]
FADE = 0.03  # fade in and out of every note [s]
PARTIALS = np.arange(1, 11)
BELL_CENTER, BELL_WIDTH = 500.0, 1.0  # the Shepard bell: center [Hz] and standard deviation [octaves]


def bell(t, f):
    """The level of a Shepard tone's partial at frequency ``f``: a bell over log frequency."""
    return np.exp(-0.5 * (np.log2(f / BELL_CENTER) / BELL_WIDTH) ** 2)


def harmonic_tone(f0, duration, fade=FADE, partials=PARTIALS):
    """An ordinary note: ten harmonics (or ``partials``) falling as 1/n, faded in and out over
    ``fade`` seconds."""
    return so.harmonic_complex(duration, FS, f0, harmonics=partials, amplitudes=1 / partials).ramp(fade)


def shepard_tone(freq, duration):
    """A Shepard tone of the pitch class of ``freq``: that pitch class in every octave from
    below 8 Hz up to Nyquist, each partial at the bell's level."""
    lowest = freq * 2.0 ** -np.ceil(np.log2(freq / 8))  # the same pitch class, between 4 and 8 Hz
    harmonics = 2 ** np.arange(int(np.floor(np.log2(FS / 2 / lowest))) + 1)
    return so.harmonic_complex(duration, FS, lowest, harmonics=harmonics, amplitudes=bell).ramp(FADE)


def sequence(notes, step):
    """Notes one after another, each ``step`` seconds after the last (they may overlap), after
    LEAD seconds of silence. ``step`` is one number, or one per note for notes of different
    lengths."""
    steps = np.broadcast_to(step, len(notes))
    onsets = LEAD + np.concatenate([[0], np.cumsum(steps[:-1])])
    length = max(int(round(onset * FS)) + note.n_samples for onset, note in zip(onsets, notes, strict=True))
    data = np.zeros(length)
    for onset, note in zip(onsets, notes, strict=True):
        first = int(round(onset * FS))
        data[first : first + note.n_samples] += note.data[:, 0]
    return so.Sound(data, FS)


def draw(snd, ax_spectrogram, ax_chroma, fmin=50, fmax=8000):
    """A spectrogram on a log frequency axis, and the sound's chroma under it."""
    stft = so.STFT(snd, win_dur=80e-3, hop_dur=10e-3)
    keep = (stft.f >= fmin) & (stft.f <= fmax)
    level = stft.db[0][keep]
    ax_spectrogram.pcolormesh(
        stft.t, stft.f[keep], level, vmin=level.max() - 70, vmax=level.max(), cmap="magma", rasterized=True
    )
    ax_spectrogram.set(yscale="log", ylabel="Frequency [Hz]", title="Spectrogram (Hann 80 ms)")
    so.chroma(snd).plot(ax_chroma)
    ax_chroma.set(xlim=(0, snd.duration))
    ax_spectrogram.set(xlim=(0, snd.duration))


def show(snd):
    """Spectrogram and chroma of one sound. Returns the figure and the panels the playhead follows."""
    fig, (ax_spectrogram, ax_chroma) = plt.subplots(2, 1, figsize=(10, 5.6), layout="constrained")
    draw(snd, ax_spectrogram, ax_chroma)
    return fig, [ax_spectrogram, ax_chroma]

bank = so.cosine_filterbank(f_lo=55, f_hi=7040, scale="cents", spacing=100)
noise = so.gaussian_noise(2.0, FS, rng=0)
back = bank.synthesize(bank.analyze(noise))
error = np.linalg.norm(back.data - noise.data) / np.linalg.norm(noise.data)
print(f"{bank.n_bands} bands, one per semitone from {bank.band_cfs[0]:.1f} to {bank.band_cfs[-1]:.0f} Hz")
print(f"Q from the formula: {1 / (2 ** (1 / 24) - 2 ** (-1 / 24)):.2f}")
print(f"analyze then synthesize, 2 s of noise: relative error {error:.1e}")

close = [so.note_to_freq(note) for note in ("C4", "E4", "G4")]
spread = [so.note_to_freq(note) for note in ("C3", "G3", "E5")]
chord = sum(harmonic_tone(f0, 2.0) for f0 in close)
power = so.chroma(chord).power[4:-4, :, 0].mean(axis=0)
order = np.argsort(power)[::-1]
names = so.views.chroma.PITCH_CLASSES
print("share of the power:", ", ".join(f"{names[i]} {power[i] / power.sum():.1%}" for i in order[:5]))

def mean_octave(snd):
    """The power-weighted mean of log2 frequency [octaves above 1 Hz], from the spectrum."""
    spectrum = np.abs(np.fft.rfft(snd.data[:, 0])) ** 2
    freqs = np.fft.rfftfreq(snd.n_samples, 1 / FS)
    keep = freqs > 0
    return np.sum(spectrum[keep] * np.log2(freqs[keep])) / np.sum(spectrum[keep])


means = [mean_octave(shepard_tone(c4 * 2 ** (step / 12), 1.0)) for step in range(12)]
print(f"mean log2 frequency over the 12 pitch classes: {min(means):.3f} to {max(means):.3f} octaves")
print(f"the bell's center, {BELL_CENTER:g} Hz: {np.log2(BELL_CENTER):.3f} octaves")

# The first two lines of Twinkle, Twinkle, Little Star, as (note, beats); each line ends on a
# note held for two beats.
TUNE = [
    ("C4", 1),
    ("C4", 1),
    ("G4", 1),
    ("G4", 1),
    ("A4", 1),
    ("A4", 1),
    ("G4", 2),
    ("F4", 1),
    ("F4", 1),
    ("E4", 1),
    ("E4", 1),
    ("D4", 1),
    ("D4", 1),
    ("C4", 2),
]
BEAT, GAP = 0.4, 0.15  # a beat, and the silence after each note [s]
melody = sequence(
    [
        harmonic_tone(so.note_to_freq(note), beats * BEAT - GAP, fade=0.06, partials=np.arange(1, 7))
        for note, beats in TUNE
    ],
    [beats * BEAT for _, beats in TUNE],
)
melody_chroma = so.chroma(melody)
played = {
    "the tune": melody,
    "its chroma, no height": melody_chroma.to_sound(FS),
    "its chroma, near E4": melody_chroma.to_sound(FS, octaves="E4"),
    "its chroma, near E6": melody_chroma.to_sound(FS, octaves="E6"),
}
for label, snd in played.items():
    measured = so.chroma(snd).power[:, :, 0]
    shares = measured.sum(axis=0) / measured.sum()
    reference = melody_chroma.power[:, :, 0].sum(axis=0) / melody_chroma.power.sum()
    difference = np.abs(shares - reference).max()
    print(f"{label:24s} largest difference in a pitch class's share of the power: {difference:.4f}")