線上教材：音樂資訊分析

時域的音高追蹤，是指直接使用 wav 檔的波形進行音高追蹤。這種方法通常會在音檔內容比較簡單的狀況下使用，例如單一樂器的演奏，或者單一人的清唱等等。本篇將由肉眼觀察波形開始，慢慢進階到半自動，乃至全自動的計算音高。部分內容在基本處理的篇章中已經介紹過，若你已經熟悉，可以自行跳過。

我們先產生一個特定頻率的正弦波，然後看看能否在假裝不知道它的頻率的情況下，只用波形來看出它的音高。以下範例會產生一個 440 Hz，長度兩秒鐘的正弦波，並且把波形圖畫出來；你可以進行局部放大，並觀察一或多個週期經過了幾個取樣點，並經由取樣頻率來推算出音訊頻率。比方說，取樣點 20,009 到 20,191（從 0 起算，下同）經過了 5 個週期，你也知道取樣頻率是 16 kHz，因此推得的音訊頻率是 16000 / ((20191 - 20009) / 5) = 439.56 Hz。雖然由於取樣等因素，而讓推算的結果跟真實的結果有些許誤差，但兩者已經相當接近。
import matplotlib.pyplot as plt
import numpy as np

DUR = 2
F = 440
FS = 16000

t = np.arange(DUR * FS) / FS
y = np.sin(2 * np.pi * F * t) * 0.8 * 32768

plt.plot(y, '.-')
plt.xlabel('Sample Index')
plt.ylabel('Raw Value')
plt.show()
如果要進行比較自動化的音高追蹤，也有許多方法，本篇先介紹一個簡單的方法，是利用自相關函數(autocorrelation function, ACF)來進行音高的自動化計算。此方法的概念是先將自己複製一份，然後不斷滑動，並將兩份之間重疊的範圍做向量內積。一開始完全重疊的時候，算出來的值會最大；而慢慢滑出去以後，內積值會漸漸變小，但是滑動距離過了半個週期，開始接近一個週期時，值又會慢慢變大，然後在一個週期的地方產生一個高峰，依此類推。由於前述步驟是針對一個 frame 來計算，因此將取樣頻率除以第一個高峰的位置，就是該 frame 的推測音高。

例如以下範例中，依然是產生了一個頻率為 440 Hz 的正弦波，但我們先框出其中一個 frame 來計算 ACF，可發現第一個高峰出現在索引值為 36 的位置，因此推測其音高為 16000 / 36 = 444.44 Hz。各位可以發現此處仍然有誤差存在，但如果仔細觀察 ACF 的計算結果，並腦補出一條曲線的話，會發現大約 36.3 的地方才是最高點，如果拿 36.3 去除 16000，計算出來就會是誤差小很多的 440.77 Hz，因此，一些比較進階的方法，會對 ACF 等方式計算出的結果進行曲線擬合，來取得誤差比較小的音高。
import matplotlib.pyplot as plt
import numpy as np

DUR = 2
F = 440
FS = 16000
FRAME_SIZE = 512
sample_start = 20000
sample_end = sample_start + FRAME_SIZE

t = np.arange(DUR * FS) / FS
y = np.sin(2 * np.pi * F * t)
y = y[sample_start:sample_end]
acf = np.array([np.sum(y[i:]*y[:FRAME_SIZE-i]) for i in range(FRAME_SIZE)])

plt.subplot(2, 1, 1)
plt.plot(y, '.-')
plt.subplot(2, 1, 2)
plt.plot(acf, '.-')
plt.show()
在上述範例中，第一高峰的位置是用目測來找尋，如果你想要自動計算的話，則不能直接取整段曲線的最大值，因為整段曲線的最大值一定會出現在開頭，而頭幾個點的位置代表過高而不合理的頻率，所以通常會將頭幾個點去掉後，再來找剩餘部分的最大值。

我們接著要來使用真人錄製的清唱聲音，來示範全自動化的 ACF 音高追蹤。以下使用的是 MIR-1K 資料集中的 geniusturtle_1_01.wav 音檔，並框出某一個特定的 frame 來進行音高追蹤：
import matplotlib.pyplot as plt
import numpy as np
import scipy.io.wavfile

WAV_PATH = 'MIR-1K/Wavfile/geniusturtle_1_01.wav'
FRAME_SIZE = 512
sample_start = 26000
sample_end = sample_start + FRAME_SIZE

fs, y = scipy.io.wavfile.read(WAV_PATH)
print(fs, y.dtype, y.shape)

y_vocal = y[:, 1]
y_vocal = y_vocal[sample_start:sample_end]
y_vocal = y_vocal / 32768
acf = np.array([np.sum(y_vocal[i:]*y_vocal[:FRAME_SIZE-i]) for i in range(FRAME_SIZE)])

max_idx = np.argmax(acf[32:]) + 32
freq = fs / max_idx
 
plt.subplot(2, 1, 1)
plt.plot(y_vocal, '.-')

plt.subplot(2, 1, 2)
plt.plot(acf, '.-')
plt.plot(max_idx, acf[max_idx], 'ro')
plt.text(max_idx+3, acf[max_idx], 'Freq: {:.4f} Hz'.format(freq))

plt.show()
在上述範例中：

我們指定的音框位置與大小，是從第 26,000 個取樣點開始，使用 512 個取樣點。

MIR-1K 資料集的音檔格式，是左聲道為背景音樂，右聲道為人聲，本範例因為只針對人聲來處理，所以只取用右聲道。

對於 ACF 計算的結果，我們扣掉前 32 個點不計，再來尋找剩下部分的最大值的位置，代表我們認為該 frame 的基礎頻率，不可能讓 ACF 的最大值出現在前 32 個點；在取樣率為 16 kHz 的情況下，相當於不會大於 16000 / (32 - 1) = 516.13 Hz。

在畫出的圖當中，上半部是該 frame 的波形，下半部是 ACF 的計算結果。

對於追蹤所有音框的音高，你可以很簡單的加上一層迴圈。我們也順便將 MIR-1K 資料集中的音高標註，跟自動追蹤的結果做對照，只是由於資料集中的標註單位是「半音」，因此需要進行適當的轉換：
import matplotlib.pyplot as plt
import numpy as np
import scipy.io.wavfile
from mir_eval.melody import evaluate
from mir_eval.util import midi_to_hz

WAV_PATH = 'MIR-1K/Wavfile/geniusturtle_1_01.wav'
LAB_PATH = 'MIR-1K/PitchLabel/geniusturtle_1_01.pv'
FRAME_SIZE = 640
HOP_SIZE = 320

with open(LAB_PATH, 'r') as fin:
	pitch_gt = list(map(float, fin.read().splitlines()))

fs, y = scipy.io.wavfile.read(WAV_PATH)
print(fs, y.dtype, y.shape, len(pitch_gt))

y_vocal = y[:, 1]
y_vocal = y_vocal / np.max(np.abs(y_vocal))

frame_time = []
frequencies = []
pitches = []
for frame_start in range(0, y_vocal.size, HOP_SIZE):
	frame_end = frame_start + FRAME_SIZE
	if frame_end > y_vocal.size:
		break
	frame_time.append((frame_start + frame_end) / 2 / fs)
	acf = np.array([
		np.sum(y_vocal[frame_start+i:frame_end+1]*y_vocal[frame_start:frame_end-i+1]) \
			for i in range(int(FRAME_SIZE/2))
	])
	max_idx = np.argmax(acf[32:]) + 32
	freq = fs / max_idx
	pitch = 69 + 12 * np.log2(freq / 440)
	frequencies.append(freq)
	pitches.append(pitch)

plt.subplot(2, 1, 1)
plt.plot(np.arange(y_vocal.size) / fs, y_vocal)

plt.subplot(2, 1, 2)
plt.plot(frame_time, pitches, '.-', label='Calculated')
plt.plot(frame_time, pitch_gt, '.-', label='Groundtruth')
plt.legend()

freq_gt = [midi_to_hz(pitch) if pitch != 0 else 0 for pitch in pitch_gt]
scores = evaluate(frame_time, freq_gt, frame_time, np.array(frequencies))
print(scores)

plt.show()
在上述範例中：

音框大小與音框間隔改為 640 與 320，是為了配合資料集的在做標註時的音框參數。

將單位為 Hz 的頻率，轉換成半音的公式，是「pitch = 69 + 12 * log2(freq / 440)」，其中的 freq 是頻率值，log2 是以 2 為底的對數，pitch 是半音值。

MIR-1K 資料集的目標是人聲的音高抽取，因此會將沒有人聲的部份的音高標記為 0。

前述的 ACF 只是時域音高追蹤的其中的一種方法，若有興趣知道其他方法，例如 NSDF、AMDF 或多種方法的綜合使用，或是以動態規劃(Dynamic Programming, DP)等方式來對結果後處理等等改進，都可以自行找資料研讀。