扒譜需要絕對音感？鋼琴自動轉錄（AMT）背後的聲音魔法

# 示範音檔：自己做一段「假鋼琴」錄音和它的標準答案 MIDI，後面每一段都用它
import numpy as np                      # 數值運算
import pretty_midi                      # 讀寫 MIDI
import soundfile as sf                  # 存 wav 檔

SR = 16000                              # 取樣率：每秒 16000 個數字

# 標準答案：(起音秒, 放開秒, MIDI 音高, 力度)
NOTES = [
    (0.0, 1.0, 48, 90),                 # C3
    (0.0, 1.0, 52, 70),                 # E3（和 C3 同時按：和弦）
    (0.0, 1.0, 55, 70),                 # G3
    (1.0, 1.5, 60, 100),                # C4（中央 C）
    (1.5, 2.0, 64, 90),                 # E4
    (2.0, 2.5, 67, 90),                 # G4
    (2.5, 3.5, 72, 110),                # C5
]


def piano_like(f0, dur, vel):
    """用「基頻 + 泛音 + 衰減」合成一個有點像鋼琴的音"""
    t = np.arange(int(dur * SR)) / SR                # 時間軸（秒）
    tone = np.zeros_like(t)
    for n in range(1, 9):                            # 第 1 到第 8 諧波（第 1 個就是基頻）
        fn = n * f0 * np.sqrt(1 + 4e-4 * n**2)       # 鋼琴的泛音會稍微偏高（不和諧性）
        if fn > SR / 2:                              # 超過取樣率一半會變成假的低頻（混疊，取樣定理）
            break
        tone += np.sin(2 * np.pi * fn * t) / n       # 越高的泛音越小聲
    env = np.exp(-3 * t)                             # 敲下去之後聲音慢慢變小
    env[-160:] *= np.linspace(1, 0, 160)             # 放開琴鍵：最後 10 ms 快速收掉
    return tone * env * (vel / 127)


def main():
    total = max(n[1] for n in NOTES) + 0.5           # 總長度多留 0.5 秒
    y = np.zeros(int(total * SR))
    for on, off, pitch, vel in NOTES:
        f0 = pretty_midi.note_number_to_hz(pitch)    # MIDI 編號 → 頻率（A4=69=440 Hz）
        s = piano_like(f0, off - on, vel)
        start = int(on * SR)
        y[start:start + len(s)] += s                 # 疊到同一條音軌上（好幾個音同時響，就是「複音」）
    y = 0.9 * y / np.abs(y).max()                    # 音量正規化，避免爆音
    sf.write("demo.wav", y, SR)

    pm = pretty_midi.PrettyMIDI()                    # 把標準答案存成 MIDI
    piano = pretty_midi.Instrument(program=0)        # program 0 = 鋼琴
    for on, off, pitch, vel in NOTES:
        piano.notes.append(pretty_midi.Note(velocity=vel, pitch=pitch, start=on, end=off))
    pm.instruments.append(piano)
    pm.write("demo.mid")
    print("已產生 demo.wav 和 demo.mid，長度", round(len(y) / SR, 2), "秒")


if __name__ == "__main__":
    main()

# 波形與取樣：聲音在電腦裡長什麼樣
import librosa                          # 音訊處理
import matplotlib.pyplot as plt         # 畫圖
import numpy as np

y, sr = librosa.load("demo.wav", sr=16000)          # 讀檔，並統一成 16 kHz
print("取樣率:", sr, "Hz")                           # 每秒幾個數字
print("總共幾個數字:", len(y), "→", len(y) / sr, "秒")

fig, ax = plt.subplots(2, 1, figsize=(10, 5))
t = np.arange(len(y)) / sr                           # 第幾個數字 → 第幾秒
ax[0].plot(t, y)                                     # 整段波形：看得出音量起伏
ax[0].set_title("Waveform")
z = slice(int(1.1 * sr), int(1.12 * sr))             # 放大 20 ms（中央 C 那段）
ax[1].plot(t[z], y[z], marker=".")                   # 每個點就是一個取樣
ax[1].set_title("Zoom 20 ms: each dot is one sample")
ax[1].set_xlabel("seconds")
plt.tight_layout()
plt.savefig("ch1_waveform.png")
print("圖存到 ch1_waveform.png")

# 頻譜圖：把聲音變成圖片（STFT、Mel、CQT）
import librosa
import librosa.display
import matplotlib.pyplot as plt
import numpy as np

y, sr = librosa.load("demo.wav", sr=16000)

# STFT：每次取 2048 點（128 ms），每 160 點（10 ms）往前移一次
S = np.abs(librosa.stft(y, n_fft=2048, hop_length=160, window="hann"))
print("STFT 形狀（頻率點 × 格數）:", S.shape)          # 1025 個頻率點
print("每個頻率點相隔:", sr / 2048, "Hz")             # 約 7.8 Hz

# Mel：把 1025 點合併成 229 個「低頻細、高頻粗」的頻帶，再取 log（dB）
M = librosa.feature.melspectrogram(y=y, sr=sr, n_fft=2048, hop_length=160,
                                   n_mels=229, fmin=30, fmax=8000)
M_db = librosa.power_to_db(M, ref=np.max)            # 能量（功率）→ 分貝
print("Mel 形狀（頻帶 × 格數）:", M_db.shape)

# CQT：每個半音一格，從鋼琴最低音 A0 開始共 88 格 → 剛好對上 88 個琴鍵
C = np.abs(librosa.cqt(y, sr=sr, hop_length=160,
                       fmin=librosa.note_to_hz("A0"), n_bins=88, bins_per_octave=12))
C_db = librosa.amplitude_to_db(C, ref=np.max)        # 振幅 → 分貝
print("CQT 形狀（琴鍵 × 格數）:", C_db.shape)

fig, ax = plt.subplots(3, 1, figsize=(10, 10))
librosa.display.specshow(librosa.amplitude_to_db(S, ref=np.max), sr=sr, hop_length=160,
                         x_axis="time", y_axis="linear", ax=ax[0])
ax[0].set_title("STFT (linear frequency)")
librosa.display.specshow(M_db, sr=sr, hop_length=160, x_axis="time", y_axis="mel",
                         fmin=30, fmax=8000, ax=ax[1])
ax[1].set_title("Log-Mel (229 bands)")
librosa.display.specshow(C_db, sr=sr, hop_length=160, x_axis="time", y_axis="cqt_note",
                         fmin=librosa.note_to_hz("A0"), bins_per_octave=12, ax=ax[2])
ax[2].set_title("CQT (one bin per semitone, 88 keys)")
plt.tight_layout()
plt.savefig("ch2_spectrogram.png")
print("圖存到 ch2_spectrogram.png")

# 泛音：鋼琴的一個音 = 基頻 + 一整排泛音
import librosa
import matplotlib.pyplot as plt
import numpy as np

y, sr = librosa.load("demo.wav", sr=16000)
S = np.abs(librosa.stft(y, n_fft=4096, hop_length=160))     # 視窗拉長，頻率看得更細
freqs = librosa.fft_frequencies(sr=sr, n_fft=4096)          # 每一列對應的頻率（Hz）
frame = librosa.time_to_frames(1.2, sr=sr, hop_length=160)  # 第 1.2 秒：只有中央 C 在響
spec = S[:, frame]

f0 = librosa.note_to_hz("C4")                               # 中央 C ≈ 261.6 Hz
print("C4 基頻:", round(f0, 1), "Hz")
for n in range(1, 6):                                       # 找第 1~5 諧波（第 1 個就是基頻）附近的最大值
    band = (freqs > n * f0 * 0.97) & (freqs < n * f0 * 1.03)
    peak = freqs[band][np.argmax(spec[band])]
    print(f"第 {n} 諧波：理論 {n * f0:.1f} Hz，實際峰值 {peak:.1f} Hz，"
          f"聽起來像 {librosa.hz_to_note(peak)}")            # 2f=C5、3f≈G5 → 八度錯誤的來源

plt.figure(figsize=(10, 4))
plt.plot(freqs, librosa.amplitude_to_db(spec, ref=np.max))
for n in range(1, 6):
    plt.axvline(n * f0, color="r", ls="--", alpha=0.5)      # 紅線：理論上的泛音位置
plt.xlim(0, 2000)
plt.xlabel("Hz")
plt.ylabel("dB")
plt.title("Middle C spectrum: evenly spaced harmonics")
plt.tight_layout()
plt.savefig("ch3_harmonics.png")
print("圖存到 ch3_harmonics.png")

# 一個音的生命週期：用示範音檔的合成方法彈一個 C4，畫出它的音量隨時間怎麼變
import librosa
import matplotlib.pyplot as plt
import numpy as np
from make_demo import SR, piano_like                 # 沿用 make_demo.py 的「假鋼琴」

HOLD = 1.2                                           # 按住 1.2 秒後放開
note = piano_like(librosa.note_to_hz("C4"), HOLD, 100)
y = np.concatenate([np.zeros(int(0.2 * SR)), note, np.zeros(int(0.5 * SR))])  # 前後留一點安靜

rms = librosa.feature.rms(y=y, frame_length=512, hop_length=160)[0]          # 每 10 ms 算一次音量
db = librosa.amplitude_to_db(rms, ref=np.max)                                 # 換成分貝
t = librosa.frames_to_time(np.arange(len(db)), sr=SR, hop_length=160)

plt.figure(figsize=(10, 4))
plt.plot(t, db)
plt.ylim(-60, 5)
on, off = 0.2, 0.2 + HOLD
plt.axvline(on, color="g", ls="--")
plt.axvline(off, color="r", ls="--")
plt.annotate("onset: hammer hits the string", (on, 0), (on + 0.08, 2), color="g")
plt.annotate("decay: slowly gets quieter", (0.8, -12), (0.55, -8))
plt.annotate("offset: key released, damper stops the string", (off, -30), (off - 0.75, -45), color="r",
             arrowprops=dict(arrowstyle="->", color="r"))
plt.xlabel("seconds")
plt.ylabel("loudness (dB)")
plt.title("Life of one piano note (C4)")
plt.tight_layout()
plt.savefig("ch3_envelope.png")
print("圖存到 ch3_envelope.png")

# NMF：把頻譜拆成「音色模板 × 什麼時候在響」
import librosa
import matplotlib.pyplot as plt
import numpy as np

y, sr = librosa.load("demo.wav", sr=16000)
C = np.abs(librosa.cqt(y, sr=sr, hop_length=160,
                       fmin=librosa.note_to_hz("A0"), n_bins=88))  # 88 鍵的 CQT
V = np.log1p(C)                                        # 壓縮大小聲差距（數值仍 ≥ 0）

# V ≈ W × H：W 是 7 個「音長什麼樣」，H 是它們「每一格多大聲」
W, H = librosa.decompose.decompose(V, n_components=7, sort=True)  # 背後是 sklearn 的 NMF；一個成分不一定剛好是一個音
print("V:", V.shape, "≈ W:", W.shape, "× H:", H.shape)

for k in range(W.shape[1]):                            # 每個成分最亮的那一格 → 對應哪個鍵
    key = int(np.argmax(W[:, k]))
    active = np.where(H[k] > 0.5 * H[k].max())[0]      # 這個成分主要在哪段時間響
    print(f"成分 {k}: 最強在 {librosa.midi_to_note(21 + key)}，"
          f"大約在 {active[0] * 160 / sr:.2f}–{active[-1] * 160 / sr:.2f} 秒")

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].imshow(W, aspect="auto", origin="lower")
ax[0].set_title("W: spectral template per component")
ax[0].set_ylabel("key (0 = A0)")
ax[1].imshow(H, aspect="auto", origin="lower")
ax[1].set_title("H: activation over time")
ax[1].set_xlabel("frame (10 ms)")
plt.tight_layout()
plt.savefig("ch4_nmf.png")
print("圖存到 ch4_nmf.png")

# Piano roll：音符表 ↔ piano roll，神經網路看的答案長什麼樣
import matplotlib.pyplot as plt
import numpy as np
import pretty_midi

FPS = 100                                             # 每秒 100 格（每格 10 ms）
pm = pretty_midi.PrettyMIDI("demo.mid")

print("音符表（起音, 放開, 音高, 力度）:")
for n in pm.instruments[0].notes:
    print(f"  {n.start:.2f}  {n.end:.2f}  {n.pitch}  {n.velocity}")

roll = pm.get_piano_roll(fs=FPS)[21:109]              # 128 個 MIDI 音高只留鋼琴的 88 鍵
roll = (roll > 0).astype(int)                         # 有在響 = 1，沒響 = 0
print("piano roll 形狀（88 鍵 × 格數）:", roll.shape)


def decode(roll, fps=FPS):
    """把 piano roll 轉回音符表：一段連續的 1 = 一個音"""
    notes = []
    for key in range(roll.shape[0]):
        row = np.concatenate([[0], roll[key], [0]])   # 頭尾補 0，方便找邊界
        diff = np.diff(row)
        starts = np.where(diff == 1)[0]               # 0→1：起音
        ends = np.where(diff == -1)[0]                # 1→0：放開
        for s, e in zip(starts, ends):
            notes.append((s / fps, e / fps, key + 21))
    return sorted(notes)


print("從 piano roll 解碼回來（起音, 放開, 音高）:")
for on, off, pitch in decode(roll):
    print(f"  {on:.2f}  {off:.2f}  {pitch}")

plt.figure(figsize=(10, 4))
plt.imshow(roll, aspect="auto", origin="lower", cmap="gray_r")
plt.xlabel("frame (10 ms)")
plt.ylabel("key (0 = A0)")
plt.title("Piano roll")
plt.tight_layout()
plt.savefig("ch5_pianoroll.png")
print("圖存到 ch5_pianoroll.png")

# 打分數：轉得好不好，用 mir_eval 怎麼算
# 拿標準答案，故意做一份「有錯的答案」，看 mir_eval 怎麼算 Precision / Recall / F1
import numpy as np
import mir_eval

# 標準答案：(起音秒, 放開秒, MIDI 音高)
ref = [(1.00, 1.50, 60), (1.50, 2.00, 64), (2.00, 2.50, 67), (2.50, 3.50, 72)]

# 模型的答案（假的，故意放四種狀況，其中三種是錯的）
est = [
    (1.03, 1.50, 60),   # 起音晚 30 ms → 還在 50 ms 內，算對
    (1.58, 2.00, 64),   # 起音晚 80 ms → 超過 50 ms，算錯
    (2.00, 2.50, 67),   # 完全正確
    (2.00, 2.50, 55),   # 多寫一個答案裡沒有的音 → 假音（FP）
]                       # C5（72）沒寫到 → 漏音（FN）


def to_mir_eval(notes):
    """把 (起音, 放開, MIDI) 轉成 mir_eval 要的格式：區間（秒）+ 音高（Hz）"""
    intervals = np.array([[on, off] for on, off, _ in notes])
    pitches = mir_eval.util.midi_to_hz(np.array([p for _, _, p in notes]))
    return intervals, pitches


ref_i, ref_p = to_mir_eval(ref)
est_i, est_p = to_mir_eval(est)

# 只看起音：音高對、起音誤差 ≤ 50 ms 就算對（offset_ratio=None 表示不看放開）
p, r, f, _ = mir_eval.transcription.precision_recall_f1_overlap(
    ref_i, ref_p, est_i, est_p, onset_tolerance=0.05, offset_ratio=None)
print(f"Note F1（只看起音）：P={p:.2f} R={r:.2f} F1={f:.2f}")
print("  手算：寫對 2 個（TP）；晚 80 ms 的 E4 和多寫的 G3 是 2 個假音（FP）")
print("  答案裡的 E4 沒被配對到、C5 沒寫，所以 2 個漏音（FN）；晚太多的音會同時算一次 FP 和一次 FN")
print("  P = 2/4 = 0.50，R = 2/4 = 0.50，F1 = 0.50")

# 動手轉一次：用現成的深度學習模型（Kong et al. 2021）把錄音轉成 MIDI
# 安裝：pip install piano_transcription_inference
# 權重 .pth 從 https://zenodo.org/record/4034264 下載，路徑當第二個參數傳進來
import sys
import librosa
from piano_transcription_inference import PianoTranscription, sample_rate  # sample_rate = 16000

audio_path = sys.argv[1] if len(sys.argv) > 1 else "demo.wav"  # 也可以換成你自己的錄音
ckpt = sys.argv[2] if len(sys.argv) > 2 else None              # 權重檔路徑（None = 自動下載）

audio, _ = librosa.load(audio_path, sr=sample_rate, mono=True) # 一樣先轉成 16 kHz 單聲道
model = PianoTranscription(device="cpu", checkpoint_path=ckpt) # 沒有 GPU 就用 cpu
result = model.transcribe(audio, "out.mid")                    # 轉錄並存成 out.mid

for n in sorted(result["est_note_events"], key=lambda n: n["onset_time"])[:10]:  # 依時間排序，印出前 10 個音
    print(f"起音 {n['onset_time']:.2f}s  放開 {n['offset_time']:.2f}s  "
          f"音高 {n['midi_note']}  力度 {n['velocity']}")
print("共轉出", len(result["est_note_events"]), "個音，存到 out.mid")

# 動手轉一次（續）：替模型轉出來的 MIDI 打分數
import sys
import numpy as np
import pretty_midi
import mir_eval

ref_path = sys.argv[1] if len(sys.argv) > 1 else "demo.mid"   # 標準答案
est_path = sys.argv[2] if len(sys.argv) > 2 else "out.mid"    # 模型的答案


def load(path):
    """MIDI → (起音放開區間, 音高 Hz)，這是 mir_eval 要的格式"""
    notes = [n for inst in pretty_midi.PrettyMIDI(path).instruments for n in inst.notes]
    intervals = np.array([[n.start, n.end] for n in notes])
    pitches = mir_eval.util.midi_to_hz(np.array([n.pitch for n in notes]))  # 注意：要用 Hz
    return intervals, pitches


ref_i, ref_p = load(ref_path)
est_i, est_p = load(est_path)

# Note F1：音高對、起音誤差 ≤ 50 ms 就算對（offset_ratio=None 表示不看放開）
p, r, f, _ = mir_eval.transcription.precision_recall_f1_overlap(
    ref_i, ref_p, est_i, est_p, onset_tolerance=0.05, offset_ratio=None)
print(f"Note F1（只看起音）：P={p:.2f} R={r:.2f} F1={f:.2f}")

# Note w/ offset：放開時間也要對（誤差 ≤ 50 ms 或音長 20%，取較大者）
p, r, f, _ = mir_eval.transcription.precision_recall_f1_overlap(
    ref_i, ref_p, est_i, est_p, onset_tolerance=0.05, offset_ratio=0.2)
print(f"Note F1（含放開）：P={p:.2f} R={r:.2f} F1={f:.2f}")