これはJulia Advent Calendar 2021の25日目の記事です。これまでにこのブログで書いたJulia Advent Calendar記事は以下の通りです。

毎年、音の分析について書いていますので、今年も音について。

楽器音のピッチ推定とYIN

数ヶ月前に、楽器音のピッチ検出（正確には基音周波数の推定）のコードを書きました。

これらの記事で紹介した方法は周波数領域で振幅ピークを見つけるというアプローチで、まあまあ正確な結果を得られましたが、計算量が多くて推定に時間がかかるのが難点でした（大きな音ファイルを一括処理しているのでしょうがないですね）。また、インハーモニシティ（倍音成分が基音周波数の整数倍からずれること）があるとうまくいかない気もします。

周波数領域ではなく時間領域で、つまりスペクトルではなく波形から基音周波数を推定する方法もあります。ピッチがある音は波形が周期性を持っているので、自己相関関数を使ってその周期を推定するというものです。単純に自己相関を使うだけだとうまくいかないときが多いので、いろいろな人が推定精度向上を目指して研究をしています。

今回は、時間領域での基音周波数推定アルゴリズムのうち、YINという方法を紹介します。YINは以下の論文で提案された手法で、自己相関関数とは少し違う方法を使っています。論文を誤読しているかもしれませんし、今回の実装に間違いがあるかもしれませんので、正しくは論文を参照してください。推定アルゴリズムの精度向上をステップごとに解説していて、丁寧かつ平易な英語なので読みやすいです。

Alain de Cheveigné and Hideki Kawahara. “YIN, a fundamental frequency estimator for speech and music.” The Journal of the Acoustical Society of America, vol.111, no.4, pp.1917–1930. April 2002. DOI: 10.1121/1.1458024

論文のSection IIに書かれた以下の6ステップに沿ってやっていきます。Section VIには機能を拡張する方法が載っていますが、今回は実装しません。

ステップ1：自己相関関数
ステップ2：差分の2乗和
ステップ3：差分二乗和関数を累積平均値で正規化
ステップ4：閾値の設定
ステップ5：二次関数による補間
ステップ6：時間方向に見回して推定精度を向上

準備

使用したバージョンは以下の通り。

julia> versioninfo()
Julia Version 1.7.0
Commit 3bf9d17731 (2021-11-30 12:12 UTC)
Platform Info:
  OS: macOS (arm64-apple-darwin21.1.0)
  CPU: Apple M1 Max
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, cyclone)

まずは、パッケージ使用のおまじないに加え、LibSndFileを使って音ファイルを読み込みます。

using FileIO: load
import LibSndFile

using Plots, LaTeXStrings
pyplot();

音ファイルを読み込んで、波形を見てみましょう。今回使うのはエレキベースの4弦の開放音（E1、約41Hz）です。

snd = load("EBass_E1.wav");
fs = snd.samplerate;
x = vec(convert.(Float64, snd.data))
t = [0:length(x)-1] / fs
plot(t, x, xlabel="Time (sec)", ylabel="Amplitude",
    ylim=[-1, +1], legend=false)

波形全体ではなく、一部分だけを抜き出して使います。窓長は100ミリ秒くらいにしました。あとで自己相関などを計算するときのために窓の外の波形も必要になるので、切り出すのは窓長の2倍の長さにします。

windowsize = nextpow(2, 0.100 * fs)   # 100ミリ秒に近い2のべき乗の長さに
startpos = 3000   # ファイル中のどこから切り出すかサンプル単位で指定
xx = x[startpos:startpos+windowsize*2-1];
txx = (0:length(xx)÷2-1) / fs
plot(txx, xx[1:length(xx)÷2], ylim=[-1, +1],
    xlabel="Time (sec)", ylabel="Amplitude", legend=false)

ステップ1：自己相関関数

信号 $\vec{x}=\left\{x _ j\right\}$ があるとき、自己相関関数は次のようになります（式1）。

$\displaystyle{ r_t\left(\tau\right) = \sum_{j=t+1}^{t+W}x_jx_{j+\tau} }$

自己相関関数が極大となる時間ズレがいくつもありますが、それらの逆数をとると倍音成分の周波数になります。そのうち周波数が最も低いものが基音周波数です。（周波数と知覚的な基音の高さが一致しない場合もありますが、ここではその議論はしません）

using StatsBase
r = autocor(xx, 0:windowsize-1);
tr = (0:length(r)-1) / fs
plot(tr, r, xlabel="Lag (sec)", ylabel="Correlation", legend=false)

ただ、高い周波数にある極大値が引っかかるといけません。例えば倍音成分とフォルマントが重なったところでは、基音よりも大きな成分になることがあります。テーパー（傾き）をつけることで、誤って高い周波数（倍音成分）を検出してしまうことを避けることができそうです（式3）。ここで $\tau _ \mathrm{max}$ （taumax）はチューニングパラメータになります（窓の大きさwindowsizeもチューニングパラメータ）。（ $\tau \leq \tau _ \mathrm{max}$ のとき。それ以外のときは $r^ {\prime\prime}=0$ ）

$\displaystyle{ r^{\prime\prime}_t\left(\tau\right) = r_t\left(\tau\right)\cdot\left(1 - \frac{\tau}{\tau_\mathrm{max}}\right) }$

# 50ミリ秒ていどのテーパー関数をつくる
taumax = nextpow(2, 0.050 * fs)
taper = ones(windowsize)
taper[1:taumax] = 1 .- (0:taumax-1)./taumax
taper[taumax+1:end] .= 0

# 自己相関関数にテーパー関数をかける
rr = r .* taper;
trr = (0:length(rr)-1) / fs
plot(trr, rr, xlabel="Lag (sec)", ylabel="Correlation", legend=false)

0.024秒のところに振幅0.65くらいのピークがあります。この逆数（1÷0.024=）41.67 Hzが基音周波数の推定値です。

以上のことは自己相関関数を使った方法の検討をしたという感じで、既往研究との比較のための記述みたいです。実はYINアルゴリズムには自己相関関数は使われていません。次からがYINのアルゴリズム。

ステップ2：差分の2乗和

自己相関関数では窓をスライドさせながら相関を計算しましたが、相関のかわりに差分の二乗和を使います。自己相関関数では相関値が高いところに基音周波数がありそうだと目星を付ける一方で、窓スライドさせながらの差分二乗和はゼロに近いほうが波形の類似性が高いということになります。そのため、以下の式を最小にする $\tau$ を探すことになります。（式6）

$\displaystyle{ d_t\left(\tau\right) = \sum_{j=1}^{W}\left(x_j - x_{j+\tau}\right)^2 }$

d = zeros(windowsize)
for tau=0:windowsize-1
    d[tau+1] = sum((xx[1:windowsize] - xx[tau+1:tau+windowsize]).^2)
end

Juliaだとfor文で回しても高速なのでうれしい。グラフを作ってみます。

td = (0:length(d)-1) / fs
plot(td, d, xlabel="Lag (sec)", ylabel=L"d_t", legend=false)

ステップ3：差分二乗和関数を累積平均値で正規化

差分が最小のところが基音周波数に関係してくるはずなのですが、このグラフの最初の値は同じ波形どうしの差分なので必ずゼロになります。「ゼロでない最小値」を探すという手もありますが、シンセサイザーのように全く同じ波形が繰り返す場合もあります。また、フォルマントの影響で高次倍音のほうに引っ張られるかもしれません。

そこで以下のように、「差分二乗和の累積平均値」で割って、低周波数に重みをつけるとともに、最初の値が最低値にならないようにします。（式8）

dd = zeros(windowsize)
dd[1] = 1
for tau=1:windowsize-1
    dd[tau+1] = d[tau+1] / (sum(d[1:tau+1])/tau)
end
tdd = (0:length(dd)-1) / fs
plot(tdd, dd, xlabel="Lag (sec)", ylabel=L"d^\prime_t", legend=false)

このグラフでディップになっているところが基音周波数の候補です。この中から最小値を探すのにfindminを使います。findminは最小値と配列の何番目にあるかを教えてくれます。

julia> findmin(dd)
(0.00439271669612706, 1154)

そのときの時間の逆数が周波数なので、以下のように基音周波数の推定ができます。

julia> 1 / tdd[findmin(dd)[2]]
41.630529054640064

もし完璧なチューニングがされていればE1音の周波数は41.2 Hzですので、あるていど推定できているといえるでしょう。

ステップ4：閾値の設定

ステップ3ではfindminを使って最小値を探しました。今回の音源ではうまく最初のディップが検出されましたが、グラフの右のほうにより深いディップがある場合も考えられます。この論文では、 $d^ \prime _ t\left(\tau\right)$ が0.1以下のディップをすべて検出し、その中で最も低い時間を基音周波数の候補としています。さらに、ディップが検出できなかった時には、単純に全体の最小値を使うとしています。

配列内の値を頭から見て行って、閾値（0.1）以下になった個所だけを抜き出せるようにします。

d_threshold = 0.1;

flagon = [];
flagoff = [];
flag = false;
for n = 1:length(dd)
    if flag == false && dd[n] <= d_threshold
        flag = true;
        append!(flagon, n);
    elseif flag == true && dd[n] >= d_threshold
        flag = false;
        append!(flagoff, n);
    end
end
if length(flagon) != length(flagoff)
    append!(flagoff, n);
end

flagonとflagoffには、閾値以下の区間の開始点と終了点が入ります。次に、それぞれの区間でのディップの最小値を計算します。

dips_at = [];
if length(flagon) == 0
    dips_at = findmin(dd)[2];
else
    for n = 1:length(flagon)
        val, ind = findmin(dd[flagon[n]:flagoff[n]]);
        append!(dips_at, flagon[n]+ind-1);
    end
end

今回の音源では結果は変わりませんが、推定された基音周波数は以下のようになります。

julia> 1 ./ tdd[dips_at[1]]
41.630529054640064

ステップ5：二次関数による補間

前のステップでは標本化された点にしか振幅がないため、計算された最小値が本当に基音周波数とぴったり一致していない場合も考えられます。そのため、最小値前後の3標本をとおる二次関数を計算して、その関数の極小値をとることにします。この計算のときは $d^ \prime$ ではなく $d$ を使うとよいとのこと。

3点を通る二次曲線の公式があるのか知らないので、愚直なやり方で解いていきます。3点 $(x _ 1, y _ 1)$ 、 $(x _ 2, y _ 2)$ 、 $(x _ 3, y _ 3)$ が与えられて、二次関数の係数 $c _ 2$ 、 $c _ 1$ 、 $c _ 0$ を求めるという、結局は連立方程式を解く課題なので、行列計算に持ち込めば簡単に解けます。

$\displaystyle{ y1 = c _ 2x _ 1^ 2 + c _ 1x _ 1 + c _ 0 \\ y2 = c _ 2x _ 2^ 2 + c _ 1x _ 2 + c _ 0 \\ y3 = c _ 2x _ 3^ 2 + c _ 1x _ 3 + c _ 0 }$

julia> n = 1;   # 試しにnが1の時だけやってみる
julia> px = [td[dips_at[n]-1] ; td[dips_at[n]] ; td[dips_at[n]+1]];
julia> py = [ d[dips_at[n]-1] ;  d[dips_at[n]] ;  d[dips_at[n]+1]];
julia> pc = py' * inv([(px.^2)' ; px' ; ones(3)'])
1×3 adjoint(::Vector{Float64}) with eltype Float64:
 3.37958e8  -1.6241e7  1.95122e5

二次関数の頂点の座標は

$\displaystyle{ \left(-\frac{c_1}{2c_2} , -\frac{c_1^2-4c_2c_0}{4c_2}\right) }$

を使って計算できます。実際のコードにまとめる前に、グラフを見てみましょう。3つのサンプル点を●で、二次補間を破線、見つけた最低点を★であらわしています。

tmpx = range(td[dips_at[n]-2], stop=td[dips_at[n]+2], length=1000)
tmpy = pc[1]*tmpx.^2 .+ pc[2]*tmpx .+ pc[3]
plot(tmpx, tmpy, label="2nd-order curve fit",
    xlabel="Lag (sec)", ylabel=L"d_t", linestyle=:dash)
scatter!(px, py, label="Points at samples", markersize=6)
scatter!([-pc[2]/(2*pc[1])], [-(pc[2]^2 - 4*pc[1]*pc[3])/(4*pc[1])],
    markersize=12, markershape=:star5, label="Minimum")

うまく補間して最小値を推定できているようです。上記ではn=1のときだけをやりましたので、すべてのnを推定するコードにします。

dips = Matrix(undef, length(dips_at), 2)
for n = 1:length(dips_at)
    px = [td[dips_at[n]-1] ; td[dips_at[n]] ; td[dips_at[n]+1]]
    py = [ d[dips_at[n]-1] ;  d[dips_at[n]] ;  d[dips_at[n]+1]]
    pc = py' * inv([(px.^2)' ; px' ; ones(3)'])

    dips[n,1] = -pc[2]/(2*pc[1])
    dips[n,2] = -(pc[2]^2 - 4*pc[1]*pc[3])/(4*pc[1])
end

その結果、ディップがある時間とそのときの $d$ の値が求まりました。

julia> dips
7×2 Matrix{Any}:
 0.0240281   2.91521
 0.0480732   8.32154
 0.072129   14.1172
 0.096186   21.1076
 0.12024    29.5267
 0.144296   39.5693
 0.168354   51.2623

これを周波数で見てみましょう。

julia> 1 ./ dips[:,1]
7-element Vector{Float64}:
 41.618007063352046
 20.801594222737307
 13.864052193957816
 10.396522087893443
  8.316702288236216
  6.930219382142619
  5.939855492937437

ちゃんと41.618 Hzと推定できていますね。

ステップ6：時間方向に見回して推定精度を向上

基音周波数の推定は瞬間値で行っていますので、音信号から取ってくる部分によっては推定を誤る場合があります。基音周波数は時間方向にあるていど定常的なはずなので、ある時刻 $t$ の前後 $T _ \mathrm{max}$ の時間範囲を見て、その範囲内での推定周波数の最低値を取り出すことになります。論文では $T _ \mathrm{max} = 25\,\mathrm{ms}$ としています（時刻$t$の前12.5 ms・後12.5 msずつです）。