Juliaで基音周波数の推定（続き）

先日、Juliaを使って基音周波数（≒ピッチ）を推定する記事を書きました。

marui.hatenablog.com

その中で、基音よりも倍音成分のほうが振幅が大きいときのために自己相関関数を使用して対処するという説明をしましたが、それでもうまくいかない音があったので、別の方法を使うことにしたよというのが今日の記事です。

手元にあったエレキベースの4弦の解放音（E1の音、約41.20 Hz）です。時間波形を見ると30秒近くまでサスティンが伸びているのがわかります。

これを先日と同じようにFFTをかけて周波数スペクトルを見てみます。

すると、第2倍音・第4倍音が大きく出ているのがわかります。以下のように前回と同じコードを使って自己相関をとって20 Hz以下を無視したグラフを描いたところ……

yc = xcorr(y, y);
yc = yc[length(yc)÷2+1 : end];
yc = yc ./ maximum(yc);

lowfreq = 20;
yc2 = yc[f .>= lowfreq];
fc2 = f[f .>= lowfreq];

plot(fc2[1:50000], yc2[1:50000],
  xlabel="Frequency (Hz)", ylabel="Correlation", legend=false)

相関値が最も高いのは基音ではなく第2倍音となってしまいました。ここでfindmax()でピーク周波数を求めると83 Hzとなり、1オクターブずれた結果になってしまいます。

自己相関関数だとうまくいかない場合もあるということで、やりかたを変えてみます。ネットで少し調べたところHarmonic Product Spectrumというテクニックがあることがわかりました。（最初の報告は A. Michael Noll. “Pitch Determination of Human Speech by the Harmonic Product Spectrum, the Harmonic Sum Spectrum and a Maximum Likelihood Estimate,” In Proceedings of the Symposium on Computer Processing in Communications, Vol. XIX, Polytechnic Press: Brooklyn, New York, (1970), pp. 779–797. らしいです。シンポジウムじたいは1969年に開催されたようです。）

下の図はリンク先からお借りしたものですが、スペクトルを周波数方向に1/2、1/3……と縮めながら積をとることで基音周波数の成分だけが突出するようになるという概念図になっています。

これをそのまま実装したのが次のコードです。DSP.resample()関数を使ってスペクトルを徐々に縮めて、ycに周波数スペクトルの積を入れていきます。第5倍音くらいまででも十分なようですが、以下では第10倍音までの積を取るようにしています。

yc = copy(y);
for n = 2:10
  yr = resample(y, 1/n);
  yr[isnan.(yr)] .= 0;
  yc[1:length(yr)] .*= yr;
end

これで計算したHarmonic Product Spectrumが以下のグラフです。

一か所だけが突出しているのがわかります。ピーク部分付近だけを拡大すると……

ここから基音周波数は41.52 Hzだと推定されました。周波数から音名を計算したところ、正確なE1からは13～14セントくらいずれていますので、チューニングをしっかりやらないとですね。

Juliaで基音周波数の推定

2021-09-24追記：正確に基音周波数が推定できない音があったので、別記事に対策を書きました。 →Juliaで基音周波数の推定（続き） - 丸井綜研

ピッチを感じる楽器の音は、基音と倍音とその他の成分音からなっています。どんな成分が入っているのか分析したいときには「楽器音の分析とスペクトログラム (Julia Advent Calendar 2018) - 丸井綜研」をぜひ。

ところで楽器音のピッチは基音の周波数と一致することが多いので、ピッチを知りたいときには基音の周波数を推定します。今日は基音の周波数を推定するためのごく簡単な方法を紹介します。

音を読み込む、処理する、可視化する、といったことをしたいので、必要なパッケージを読み込みます。

using FileIO: load
import LibSndFile
using DSP, FFTW
using Plots
pyplot();

音信号をファイルから読み込むのにはLibSndFile.load()を使います。今回はクラシックギターの1弦の開放音を読み込みます。

snd = load("guitar.wav")

読み込んだ信号の波形を見てみましょう。波形はsnd.data、標本化周波数はsnd.samplerateに入っていますが、名前が長いのでそれぞれxとfsに代入して使うことにします（モノラルファイルではない場合には1チャンネル目だけを使います）。tは各サンプルの時刻を秒単位であらわしたものです。

x = snd.data;
if length(size(x)) > 1
  x = x[:,1];
end
fs = snd.samplerate;
t = range(0, length(x)-1, step=1) ./ fs;
plot(t, x, ylim=[-1, +1],
  xlabel="Time (sec)", ylabel="Amplitude", legend=false)

一回の撥弦音が収録されているのがわかります。

倍音成分を見るために、フーリエ変換を使って周波数スペクトルの計算をします。普通のFFTを使ってもいいのですが、負の周波数成分が出てこなくて手軽なFFTW.rfft()を使い、振幅を計算します。正確な振幅値を得るためには信号長で補正をかけたりするんですが、今回は振幅の絶対的な値には興味がないので無視します。fはスペクトルの各周波数ビンの中心周波数です。（普通のFFTを使った分析については「パワースペクトル計算の2つの方法 (Julia Advent Calendar 2019) - 丸井綜研」に書きました）

y = abs.(rfft(x));
f = range(0, stop=fs/2, length=length(y));
plot(f, y,
  xlabel="Frequency (Hz)", ylabel="Amplitude", legend=false)

低い周波数に何本かピークがあるのがわかります。拡大表示してみましょう。

plot(f[1:50000], y[1:50000],
  xlabel="Frequency (Hz)", ylabel="Amplitude", legend=false)

倍音成分が330 Hz付近、660Hz付近、990 Hz付近……に出ています。基音の成分が突出しているので、振幅が最大値になる周波数を求めれば基音周波数が推定できそうです。

(val, ind) = findmax(y)
f[ind]

とすると、330.96 Hzという結果でした。A4を440 Hzに合わせたとするとE4は329.6 Hzですから、だいたい正しい周波数が推定できてそうです。しかし、楽器によっては基音よりも第2倍音や第3倍音が振幅が大きい場合もあります。その場合、上記のようにスペクトルの最大値だけを見てしまうと、基音の周波数ではなく第2倍音や第3倍音の周波数を検出してしまいます。

そこで、スペクトルの自己相関関数を計算します。ピッチを持つ楽器では多くの場合、倍音成分は基音の整数倍の周波数にあらわれます。周波数軸上で等間隔に並んでいるので、自己相関を計算すればその周期を検出できるわけです。自己相関の計算には相互相関関数を使いますが、結果が対象になるので、片側半分だけを取り出します。

yc = xcorr(y, y);
yc = yc[length(yc)÷2+1 : end];
yc = yc ./ maximum(yc);
plot(f[1:50000], yc[1:50000],
  xlabel="Frequency (Hz)", ylabel="Correlation", legend=false)

ちょうど330 Hz付近にピークがあります。自己相関なので当然ながらシフト量がゼロのときが相関最大になります。これを除去しないとfindmax()が使えませんので、下記のようにlowfreqという変数を用意して、その周波数以上の成分だけを検出対象にします。

lowfreq = 20;
yc2 = yc[f .>= lowfreq];
fc2 = f[f .>= lowfreq];
(val, ind) = findmax(yc2);
plot(fc2[ind-1000:ind+1000], yc2[ind-1000:ind+1000],
  xlabel="Frequency (Hz)", ylabel="Amplitude", legend=false)

上のグラフは相関最大になった付近の拡大図です。もとのファイルの長さがけっこうあったので、FFTで得られたビン数も多く周波数分解能も高くなっています。ピークの周波数を見てみます。

fc2[ind]

結果、330.28 Hzと、基音周波数の推定ができました。

Juliaで手軽に基音周波数の推定をする簡単な方法を紹介しました。自分が手元で使っているコードでは、入力信号が短い場合はゼロ詰めして周波数分解能を上げるという処理を加えたものを関数化しています。

基音周波数推定は奥が深く、リアルタイム動作に適した洗練された方法も提案されていますので、興味がある方は調べてみてください。（たとえばDOI: 10.1121/1.1458024とかが有名です。「Juliaで基音周波数の推定（YINの解説）（Julia Advent Calendar 2021） - 丸井綜研」で解説してみました）

本日の給油

給油日	オドメーター (km)	給油量 (L)	単価 (円/L)	燃費 (km/L)	距離単価 (円/km)
2020-09-21	11760.8	3.04	130.59	48.88	2.67
2021-02-21	11900.9	2.92	139.73	47.98	2.91
2021-03-12	12042.2	3.08	133.77	45.88	2.92
2021-05-03	12177.4	2.46	139.84	54.96	2.54
2021-06-06	12290.5	1.99	142.21	56.83	2.50
2021-07-31	12417.3	2.62	148.09	48.40	3.06

ここ1ヶ月のうちにバイクの点検をしてもらった。ペダルのゴムみたいな消耗品の交換に加えて、オイル交換してタイヤの空気圧を調整するていど。そろそろ9年になるけど丁寧に乗っているね、とのこと。

今日は、近所のバイク販売店を回って新車を見たり（バイクレンタル会員に申し込んでみた）、夕食の食材を仕入れてバターチキンカレーを作ったなど。

ホットクックを使ってスパイスカレーを作った

はじめに

料理そのものは嫌いじゃなくて、時間があれば気楽にパパっと作りたいとは思ってます。いっぽうで時間をかけてゆっくりと作るものも好きです。たとえば発酵食品作りみたいなのも好きで、イースト菌や麴を使ってごにょごにょやってみたり、ヨーグルトを作ったりしてました。発酵は基本的に放置でいけますし、多くても一日一回確認すればいいかんじなので気楽です。

一方で、手を出してはいけない最後の聖域は、蕎麦打ち、包丁研ぎ、スパイスカレーだと思っていました。無限に探究できてしまいそうで、かなりヤバい。（製麵機を買いそうになったときもありました）

www.seimen.club

とにかく手軽なスパイスカレー

とくにカレーは無数にあるスパイスの沼にはまるのが目に見えています。ところが、どこかで印度カリー子さんによる「クミン、ターメリック、コリアンダーだけあれば、スパイスカレーが始められる」という記事を読んだんですね。それが頭のどこかにあって、Kindleで読む本をAmazonで見ていたら、印度カリー子「ひとりぶんのスパイスカレー」（https://amzn.to/2ToNKPy）がPrime Readingで無料だったんです。読むだけだから～とポチって読み始めたらやっぱり面白そう。

www.hotpepper.jp

カレーのベースになるグレイビーを作っておけば、あとは食べたいときに手軽にカレーが作れるということで、6月上旬にカレーを2回ほど作りました。とてもおいしくできたんですが、まず作らないといけないグレイビーのタマネギ炒めが結構大変で、その部分をなんとかしたくなりました。

ホットクック

うちにはホットクック（https://amzn.to/3h9M31k）があるので、それを使って大量にグレイビーを作り置きしておけば、いつでもカレーが食べられます。似たようなことを考えている人は結構いるみたいで、時短スパイスカレー作りにホットクックを使っている人を見つけました。

note.com

で、その記事から参照されている、ホットクックで飴色タマネギを作る方法が以下のサイト。

kitchen-panda.com

記事で使っていたホットクックは大型のほうなので、うちの小型のほうでは同じ分量と時間ではできないかもしれません。とりあえず「炒める」「煮詰める」という方法でいけることがわかったので、以下に手順をメモしておきます。

本日の給油

給油日	オドメーター (km)	給油量 (L)	単価 (円/L)	燃費 (km/L)	距離単価 (円/km)
2020-03-16	11037.1	3.46	138.73	42.51	3.26
2020-03-20	11217.6	2.83	133.92	63.78	2.10
2020-03-21	11407.0	3.11	142.12	60.90	2.33
2020-03-22	11612.2	3.50	138.57	58.63	2.36
2020-09-21	11760.8	3.04	130.59	48.88	2.67
2021-02-21	11900.9	2.92	139.73	47.98	2.91
2021-03-12	12042.2	3.08	133.77	45.88	2.92
2021-05-03	12177.4	2.46	139.84	54.96	2.54
2021-06-06	12290.5	1.99	142.21	56.83	2.50

キャンプ用品を見に行ったり、バイクを見に行ったり、夕食の食材を仕入れたり、バイク洗車したりした週末（ピカピカの展示車を見て、汚れたままのバイクでバイク屋さんに行ったのが恥ずかしくなってしまったのでした）。ガソリン単価が上がってきてますね。

【追記】6月27日(日)に総走行距離が12345.6 kmになりました。

本日の給油

給油日	オドメーター (km)	給油量 (L)	単価 (円/L)	燃費 (km/L)	距離単価 (円/km)
2020-03-16	11037.1	3.46	138.73	42.51	3.26
2020-03-20	11217.6	2.83	133.92	63.78	2.10
2020-03-21	11407.0	3.11	142.12	60.90	2.33
2020-03-22	11612.2	3.50	138.57	58.63	2.36
2020-09-21	11760.8	3.04	130.59	48.88	2.67
2021-02-21	11900.9	2.92	139.73	47.98	2.91
2021-03-12	12042.2	3.08	133.77	45.88	2.92
2021-05-03	12177.4	2.46	139.84	54.96	2.54

ゴールデンウィーク中にいろいろと旅行したかったんだけど、COVID-19状況下ではそうも行かず。職場からは「飲食をともなう会合は禁止」「県境をまたぐ移動は自粛」という指示があったので、県境をまたがないソロ移動ということで、羽田の福太郎で明太子を買い、明治神宮で交通安全御守をいただいてきた。

youtu.be

AmbisonicsとUnityの座標系

Ambisonicsの本を読んでいて、頭が悪いので座標系がわけわからなくなりがちだったので、絵に描きました。Ambisonicsは前後がX軸・左右がY軸・上下がZ軸。コンピュータ画面の座標系（左上が原点で右がX軸・下がY軸）や、Unity（右左がX軸・上下がY軸・前後がZ軸）と違うので混乱していたのでした。右手系・左手系の違いに加えて、上下軸の名前が変わるのもややこしい。*1

読んでいたAmbisonicsの本はこれ。MS方式から入ってHigher Order Ambisonicsまで説明があるんですが、すべての章においてフリーソフトでどうやって実現するのかという話題が書かれていて、とても実践的です。

www.springer.com

平易な英語で書かれているし、電子版は無料でダウンロードできます。Amazonには無料のKindle版もあるので、興味のある人にはおすすめです。（Amazon.co.jp）

*1:さらにややこしいことに、Blenderは（起動画面では）後前がX軸・右左がY軸・上下がZ軸の右手系で、Unreal Engineは（サードパーソンプロジェクトの起動画面では）前後がX軸・右左がY軸・上下がZ軸の左手系、（2Dサイドスクロールプロジェクトでは）右左がX軸・上下がZ軸となっています。Unrealは主人公が進む方向がX軸になっているのかも。