Keyboard Piano

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132

#!/usr/bin/env python3
"""Regenerate samples/<midi>.ogg for Keyboard Piano.

A physical-ish piano model in pure Python, encoded to Ogg Vorbis by ffmpeg:

- every note is a stack of stiff-string partials, stretched sharp by an
  inharmonicity that grows up the keyboard, as many as fit under Nyquist
- the hammer strikes an eighth of the way along the string, which notches
  every eighth partial, and its felt rolls the spectrum off
- each partial decays in two stages, a fast attack decay and a long
  aftersound, and upper partials die sooner
- two or three strings per note, unequal and slightly out of tune, so the
  tone shimmers without ever beating down to silence
- a few milliseconds of filtered noise for the hammer knock

48 kHz mono, so the usual 48 kHz output device plays it without resampling.

    python3 tools/synth.py            # writes ./samples, needs ffmpeg
    python3 tools/synth.py out_dir
"""

import array
import math
import multiprocessing
import os
import random
import subprocess
import sys
import tempfile
import wave

RATE = 48000
LOW, HIGH = 36, 91  # C2..G6: the key map's range across every octave shift
HAMMER_AT = 0.125
PEAK = 0.85
VORBIS_QUALITY = "4"


def note_seconds(midi):
    # bass strings ring longest: 4.0 s at the bottom, 1.5 s at the top
    return 4.0 - 2.5 * (midi - LOW) / (HIGH - LOW)


def add_partial(out, freq, amp, t60, phase):
    """Add amp * exp(-t/tau) * sin(w t + phase) by a two-multiply recurrence."""
    w = 2 * math.pi * freq / RATE
    r = math.exp(-math.log(1000) / (t60 * RATE))
    c, r2 = 2 * r * math.cos(w), r * r
    y2 = amp * math.sin(phase - w) / r
    y1 = amp * math.sin(phase)
    for i in range(len(out)):
        out[i] += y1
        y1, y2 = c * y1 - r2 * y2, y1


def render(midi):
    rng = random.Random(midi)
    f0 = 440.0 * 2 ** ((midi - 69) / 12)
    n = int(RATE * note_seconds(midi))
    out = [0.0] * n

    up = (midi - LOW) / (HIGH - LOW)  # 0 at the bottom of the range, 1 at the top
    stiffness = 2e-4 * 2 ** ((midi - 48) / 12 * 1.2)
    felt = 2400 + 1800 * up  # Hz: where the hammer's spectrum has fallen by 1/e
    sustain = 11.0 * 2 ** (-2.6 * up)  # seconds for the fundamental's aftersound to fall 60 dB
    strings = [(0.0, 1.0)] if midi < 40 else [(0.0, 1.0), (0.7, 0.7)]
    if midi >= 48:
        strings.append((-1.0, 0.5))

    k = 1
    while k <= 64:
        stretched = k * f0 * math.sqrt(1 + stiffness * k * k)
        if stretched > RATE * 0.42:
            break
        amp = abs(math.sin(k * math.pi * HAMMER_AT)) * math.exp(-stretched / felt) / k**0.8
        if amp > 1e-4:
            t60 = sustain / (1 + 0.14 * (k - 1))
            for cents, weight in strings:
                freq = stretched * 2 ** (cents / 1200)
                phase = rng.uniform(0, 2 * math.pi)
                add_partial(out, freq, amp * weight * 0.7, t60 * 0.3, phase)
                add_partial(out, freq, amp * weight * 0.3, t60, phase)
        k += 1

    # the hammer knock: a short burst of low-passed noise
    body = max(abs(s) for s in out)
    smooth = math.exp(-2 * math.pi * (900 + 2500 * up) / RATE)
    low = 0.0
    for i in range(int(RATE * 0.03)):
        low = smooth * low + (1 - smooth) * rng.uniform(-1, 1)
        out[i] += 0.5 * body * low * math.exp(-i / (RATE * 0.006))

    attack, release = int(RATE * 0.002), int(RATE * 0.08)
    for i in range(attack):
        out[i] *= 0.5 - 0.5 * math.cos(math.pi * i / attack)
    for i in range(release):
        out[n - 1 - i] *= i / release

    scale = PEAK * 32767 / max(abs(s) for s in out)
    return array.array("h", (int(s * scale) for s in out))


def write_note(job):
    midi, out_dir, tmp = job
    pcm = render(midi)
    if sys.byteorder == "big":
        pcm.byteswap()
    wav = os.path.join(tmp, f"{midi}.wav")
    with wave.open(wav, "wb") as f:
        f.setnchannels(1)
        f.setsampwidth(2)
        f.setframerate(RATE)
        f.writeframes(pcm.tobytes())
    subprocess.run(
        ["ffmpeg", "-v", "error", "-y", "-i", wav, "-c:a", "libvorbis", "-q:a", VORBIS_QUALITY,
         "-map_metadata", "-1", "-fflags", "+bitexact", os.path.join(out_dir, f"{midi}.ogg")],
        check=True,
    )
    return f"{midi}.ogg  {len(pcm) / RATE:.2f} s"


def main():
    out_dir = sys.argv[1] if len(sys.argv) > 1 else "samples"
    os.makedirs(out_dir, exist_ok=True)
    with tempfile.TemporaryDirectory() as tmp, multiprocessing.Pool() as pool:
        jobs = [(midi, out_dir, tmp) for midi in range(HIGH, LOW - 1, -1)]
        for line in pool.imap_unordered(write_note, jobs):
            print(line)


if __name__ == "__main__":
    main()