頭部伝達関数 (HRTF: Head-Related Transfer Function) を使って8D立体音響をPython で実装してみました。8D立体音響というのは正直言って私も厳密な定義はわかりませんが、この記事で言っているのはYoutube によくある頭の周りをクルクルまわる音源のことです。「左右のチャネルの音量を変えれば作れるのでは?」という声を聞きますが、それではできず、頭部伝達関数 というものを使う必要があります。
注意点として、8D立体音響を作るソフトがどのように処理しているかはわからないため、細かいところは異なることをご承知おきください。
頭部伝達関数 とはヒトが音の方向を知覚するために重要な役割を果たすものです。
参考文献 [1] では、頭部伝達関数 を以下のように定義しています。
音波は鼓膜に届く直前に頭や耳介,あるいは胴体の影響を受ける。このような,頭部周辺による入射音波の物理特性の変化を周波数領域で表現したものを頭部伝達関数 という。
つまり、頭や耳、鼻などが与える音への影響のことです。
例えば、図1のように左側に音源があるとき、音源と左耳の間には何もないため、音はあまり変化なく左耳に伝わります。一方、音源と右耳の間には頭や鼻という障害物があるため、右耳に伝わる音は減衰したり、左耳と比べると遅れて音が届いたりします。
図1:音が伝わる模式図
このように左右の耳で異なる音の伝わり方をすることで、ヒトは音の位置や方向を知覚できるそうです。冷静に考えると驚くハナシで、信号処理では音が到来する2次元方向を把握するのも結構難しいので、音源の3次元の方向や位置がわかる人間はすごいなと思います。
8D立体音響を作成する場合は、測定した頭部伝達関数 を音データに作用させることで左右前後に音源があると錯覚させることができます。
頭の形などは個人差があるので、他人の頭部伝達関数 を使うと前後方向については上手く錯覚させることはできませんが、左右方向については上手く錯覚させることができます。
頭の周りをまわる音源の作成
Youtube によくある頭の周りをまわる音源を作成したいと思います。
音源の位置の計算
音源が頭の周りをまわっている様子を図2に示します。
図2:音源の位置の計算
音源が図2のように角速度 [deg/s] でまわっているとき、サンプル番号 n のときの音源方向角度 [deg] は以下のように求まります。
ここで、 はサンプリング周波数です。
音源方向角度が360度以上の場合、プログラム上扱いにくいので、0度以上360度未満になるように360度の倍数で引きます。
今回使用するHRTFのデータベースは5度間隔でHRTFが測定されています。そのため、音源方向角度が7度や8度のような中途半端な場合、どのようなHRTFを使えばいいかという問題があります。
今回は、参考文献 [2] に記述されている線形2点補間 を用いて、HRTFの補間を行います。計算は簡単で以下のようにHRTFを補間します。
ここで、 は の2つのHRTFの内分比を表します。例えば、音源方向角度が8度の場合、 は5度のHRTF、 は10度のHRTFを用いて、 は0.4とします。
補足 :今回行う線形2点補間は参考文献 [2] のやり方とは少し異なっています。参考文献 [2] ではHRTFの振幅応答を補間して求めて、位相応答はHRTFが最小位相フィルタとなるように求めています。私の補間の仕方では位相が打ち消しあって、変な音になる可能性がありますが、処理した音はそこまで問題がなかったのでこのままにしています。
オーバーラップ加算法
いままでデータベースの中にはHRTFが入っているような言い方をしていましたが、厳密には頭部インパルス応答(HRIR:Head-Related Impulse Response) が入っています。HRIRをFFT することで、HRTFになります。
HRIRを波形データに畳み込めば、クルクルまわる音源ができますが、かなり時間がかかってしまいます。また、for文の処理が遅いPython ではなおさら計算時間がかかります。
そこで、参考文献 [1] に記述されているオーバーラップ加算法 を行います。オーバーラップ加算法の手順は以下です。
(1) 音源信 号 x からchg_len ずつずらしながらHRIRの長さ512点を取り出していく。
(2) 取り出した音源信 号とHRIRに512点の0埋めを行う
(3) FFT 、複素乗算、逆FFT を施してデータ長 1024点の結果を得る。
(4) 処理した結果を加算していって、出力信号 y とする。
図3:オーバーラップ加算法
プログラム
8D立体音響を作成するプログラム 3Dsound.py は以下です。
import soundfile as sf
import numpy as np
from scipy.fft import rfft, irfft
import scipy.signal as sg
wav_name = "sound_mono.wav"
out_name = "soundout.wav"
elev = 0
hrtf_dir = "hrtfs/elev" +str (elev)+"/"
N = 512
chg_len = 128
omega = 30
x, fs = sf.read(wav_name)
if len (x.shape) == 2 :
x = x[:,0 ]
x_len = len (x)
x = np.pad(x,[0 ,N*2 ], "constant" )
hrir_l = np.zeros((72 , N))
hrir_r = np.zeros((72 , N))
for i, angle in enumerate (range (0 ,360 ,5 )):
angle_str = str (angle)
angle_str = angle_str.zfill(3 )
path = hrtf_dir + "L" +str (elev)+"e" + angle_str + "a.dat"
with open (path) as f:
for j, s_line in enumerate (f):
hrir_l[i,j] = float (s_line)
path = hrtf_dir + "R" +str (elev)+"e" + angle_str + "a.dat"
with open (path) as f:
for j, s_line in enumerate (f):
hrir_r[i,j] = float (s_line)
HRTF_L = np.zeros((72 , N+1 ), dtype=np.complex128)
HRTF_R = np.zeros((72 , N+1 ), dtype=np.complex128)
for m in range (72 ):
h = np.pad(hrir_l[m,:], [0 ,N], 'constant' )
HRTF_L[m,:] = rfft(h)
h = np.pad(hrir_r[m,:], [0 ,N], 'constant' )
HRTF_R[m,:] = rfft(h)
y = np.zeros((len (x),2 ))
n_frame = x_len // chg_len + 1
for i in range (n_frame):
theta = omega * i * chg_len / fs
while int (theta) > 359 :
theta = theta - 360
m = theta / 5
m1 = int (m)
m2 = m1 + 1
if m2 == 72 :
m2 = 0
r2 = m - int (m)
r1 = 1.0 - r2
x_N = np.pad(x[i*chg_len:i*chg_len+N], [0 ,N], 'constant' )
X = rfft(x_N)
YL = X * (r1*HRTF_L[m1,:]+r2*HRTF_L[m2,:])
YR = X * (r1*HRTF_R[m1,:]+r2*HRTF_R[m2,:])
y[i*chg_len:i*chg_len+2 *N, 0 ] += irfft(YL)
y[i*chg_len:i*chg_len+2 *N, 1 ] += irfft(YR)
y = y/np.max(y)
sf.write(out_name, y[:x_len,:], fs, subtype="PCM_16" )
6~13行目 :使用するHRTFの仰角やchg_len、角速度などのパラメータを入力しています。
15~20行目 :WAVファイルから波形データを読み込んでいます。WAVファイルがステレオデータの場合は、左チャネルだけ取り出します。
22~37行目 :0度から355度のHRIRをnumpy配列に読み込んでいます。
39~46行目 :読み込んだHRIRを0埋めして、FFT することでHRTFを作成しています。
52~55行目 :フレームごとに移動音源の角度を求めています。
56~63行目 :HRTFを線形2点補間するためのパラメータを求めています。例えば、音源方向角度が17度の場合、17/5=3.4 でm1=3, m2=4 となります。インデックス3と4に対応するHRIRは15度と20度です。また、r2=3.4-3=0.4 、r1=1-r2=0.6 となり、線形補間のパラメータが求まります。
64~66行目 :x から512点取り出して、0埋めをして、FFT をしています。
67~72行目 :補間して求めたHRTFとXを乗算して、逆FFT をして、加算しています。
74~76行目 :ノーマライズ して、波形データをファイルに書き込んでいます。
処理結果
頭部伝達関数 を使って作成した8D立体音響は以下のようになります。仰角は0度、chg_lenは128点で作成しました。ヘッドホンやイヤホンでお聴きください。
入力音源
処理結果(角速度ω:30 [deg/s])
処理結果(角速度ω:120 [deg/s])
そんなに雑音もなく作成できていると思います。角速度ω:120 [deg/s] はさすがに少し気持ち悪くなりますね。
おわりに
HRTFを使って8D立体音響を実装してみました。Youtube の8D立体音響ではボーカルだけこのような処理をして、他の伴奏とかは普通にミックスしているのかもしれません。あとはコンサートホールなどのインパルス応答も畳み込んだりして、いい感じにしているのだろうと思います。
参考文献
[1] 飯田 一博、”頭部伝達関数 の基礎と3次元音響システムへの応用”、コロナ社 、2017.
[2] 西野隆典, 梶田将司, 武田一哉, 板倉文 忠, "水平方向及び仰角方向に関 する頭部伝達関数 の補間," 日本音響学会 誌, 57巻, 11号, pp.685-692, 2001.
使用したデータについて
この記事で信号処理した楽曲は Cambridge Music Technology で提供されている The Balazs Daniel Boogie Woogie Trio の Own Way To Boogie を使用させていただきました。
【楽曲を提供している Cambridge Music Technology のページ】
https://www.cambridge-mt.com/ms/mtk
頭部伝達関数 につきましては西野隆典研究室のホームページにある頭部伝達関数 データベースからHRTF data (2)を使用させていただきました。
【頭部伝達関数 データベースのページ】
http://www.sp.m.is.nagoya-u.ac.jp/HRTF