MATLABでHRTF～頭部伝達関数とは？～聴覚の仕組みは解明されていない～
R2023bで SOFA ファイルを直接取り扱えるようになり、外部 API を使用する必要がなくなりました。
その件を記事中に追記しています。
2023/09/18 更新
電子書籍発売しました！（2023/01/22）
動画・音声リンク付き紹介記事はこちら。
Amazon Kindle 電子書籍「MATLAB で簡単オーディオ プラグイン開発」
スクリプトはどなたでも無料でダウンロードできます。
MATLAB による VSTプラグイン開発の記事は以下リンクから
・MATLABでVSTプラグイン開発 [初級編]～数分？で作れるVST～LMSで無相関プラグイン
・夏休みはVST！ MATLABで楽々（？）VSTプラグイン開発
HRTFを使わない、イヤホン・ヘッドホン・ネックバンドスピーカー向け頭外化処理技術 AudiiSion EP の最新デモ音源はこちら
HRTFを使わない、狭間隔スピーカー向け空間音響再現技術 DnoteLR+の開発秘話はこちら
2022/07/22 更新
前書き
さて最近、今まで一般の場ではほとんど聞かれることのなかった、「HRTF（頭部伝達関数）」という単語がちらほらツイッター等でも見かけるようになりました。
そう、あのゲーム機の影響です。
私はHRTFが専門というわけではなく、自分のHRTFを測定したこともありませんが、何度か試してみたことはあるので、簡単な使い方を解説してみようかと思います。
HRTFとは
HRTF自体は1970年代前半からある古い技術です。
音は音源の位置によって、左右の耳に届く時間や音量がわずかに違ったり、胴体で反射したり、耳介のくぼみ等で共鳴したり小さく反射したり、横方向からであれば、逆側の耳へは頭で遮られたり、頭を回り込んだり、様々な経路を通って、その人の体の影響を受けてから耳に届きます。
人間は（他の動物も）、これらの影響の違いを学習して、音源方向を推定していると言われています。
「これらの影響」を音源から耳への経路の特性として表したのが、頭部伝達関数（Head Related Transfer Function）です。
この特性はかなり複雑で、体の影響を受けているので個人個人全く違っているという特徴があります。
ですので、もし他の人のHRTFが掛かった音を聞くと、音源方向を正しく認識できなかったり、音質が大きく変わってしまったりします。
また、耳介を粘土等で形状を変えると、最初は音源方向が曖昧になりますが一ヶ月もすれば正しい定位認識ができるようになり、元に戻しても定位認識はできる、と言う報告※1もあります。その論文では、「上書きではなく、新しい言語を覚えるように両方のHRTFを認識できるようになる」とあります。
※1 Paul M. Hofm, Jos G.A. Van Riswick, John Van Opstal: Relearning sound localization with new ears. Nature Neuroscience, 1998, 1(5), 417-421.
HRTFの周波数特性の違いの例を、以下にグラフで示します。
音源位置（水平角度）による違い
人による違い
同一人物の左耳、右耳による違い
実はステレオスピーカーで音を聞く場合は、「ちゃんと」セッティングしてあげれば、普通のステレオ音源でもかなり立体的な音像として聴くことができます。
これは、「自分のHRTF」が掛かって聴いていることも大きく影響していると考えられます。
イヤホンやヘッドフォンの場合、ほぼHRTFが掛からない状態の音が耳に入るため正常な音像定位認識が行えず、左右の耳の間の狭い範囲に音像が押し込められてしまいます（頭内定位）。
（ちなみに専門用語的には、ヘッドバンド付きのイヤホンをヘッドフォンと呼ぶので、イヤホンが上位概念です）
仮想音像再現
そこで考えられるのが、擬似的にHRTF処理を掛けた音をイヤホンから流す方法です。
手っ取り早いのは、左右の耳にマイクを付けて録音し、その位置で音を再生することです。
これがバイノーラル録音と言われるものですが、なにしろHRTFは人によって違い、少しでも合わないと効果を感じられなかったり、特に前後を誤認識したり、音質が大きく変わってしまったりします。
つまり、正確に音場を再現するには、事実上本人を使って録音しなければならないことになります。
バイノーラルの歴史は古く、1881年パリ電気博覧会が初公開と言われており、実はいわゆる「ステレオ」システムより半世紀も前です。
もう一つは、音源をマルチチャネルあるいはオブジェクト毎に録音しておいて、再生時に、聞く人のHRTFを用いて処理を掛ける方法です。
ただ、HRTFを正確に測定するのはかなり手間が掛かり、オブジェクトごとに違うHRTFを掛ける必要があるので処理も重くなります。
最近は耳の写真をスマホで撮って個人のHRTFを作成できるシステムも出てきましたが、片耳を横から撮った写真一枚から推定するようなので、簡易的なものなのかなとは思います。
HRTFは左右でかなり異なりますし、正面から見た耳介の角度、胴体の反射等も影響するはずですから。
ちなみに、「私、定位認識が不得手なんですよね」という方、鏡で自分の耳介の開き具合をチェックしてみてください。左右どちらか、あるい両方の耳介が寝てませんか？
その場合はそれを起こしてあげると、定位が認識しやすくなると思います。
それか付け耳を。
私は右耳が寝てるせいか、イヤホンでも右の定位認識がイマイチです。
少し話が逸れましたがまとめると、HRTFは音源から耳までの物理特性であり、HRTFが正確に再現されれば 360度音場が再現可能とされています。しかしHRTFは個人差が大きく、特性が複雑で処理に必要な演算量も多く、測定自体も難しく、正確に再現することは難しいのです。
そして、正確に再現されない限り、その後の脳の認識は大きく違ってしまいます。
聴覚の仕組みは解明されていない
実は、聴覚の仕組みはまだよく分かっていません。
HRTF自体もどの部分がどう影響しているのか正確には分かっていませんし、聴覚はさらに複雑です。
音像認識
ドラマや映画を見るとき、スピーカーで聴いていてもイヤホン／ヘッドフォンで聴いていても、俳優の台詞はちゃんとその口元から聞こえてきませんか？
TVでチャイムが鳴ると、つい玄関の方を見てしまいませんか？
同じところから音が出ていても、鳥の声や雷は上から聞こえたり、地鳴りやコインが床に落ちる音は下から聞こえてきたり、ささやき声は近くに、叫び声は遠くに聞こえることなども知られています。
つまり音像定位認識は、視覚や経験の影響を受ける、脳のかなり高次の認識であり、同じところから音が出ていても、それを音の中身によって空間に再配置する音像空間認識能力が人間には備わっているということです。
そうなると、猿と人間でさえ仕組みが違う可能性も大きいので、なかなか解明は難しいのでしょう。
最近の生理医学系の論文でも、
「新しい聴覚・定位認識モデルを提案した」
→ 「さらなる研究が必要である」
というのがいくつも見られます。
みんな違って聞こえる
仮想音像空間再現で一番問題なのは、どんなに耳の良い人が立体的に聞こえるように作り込んでも、他の人には全く同じようには聞こえないということです。
音が他の人にどう聞こえているのかはなかなか分かりません。
映像と違って、「ほらここが」と確認することもできません。
私は「標準」のHRTFとかなり違うのか、今まで聴いてきたバーチャルサラウンド・頭外定位系は、一度も謳い文句通りに聞こえたことがありません・・。
これはなにも開発者が嘘を言っているわけではなく、本人にはそう聞こえているのでしょう。
今まで学会等でのデモを通じて数百人の方々と色々お話しさせて頂きましたが、「謳い文句通りに聞こえない」人はかなり多いという印象です。
開発側の、「こう聞こえて欲しい！」という希望的観測も入っているのかも・・。
それと、聴く側にも「問題」があるのかもしれません。
音像に限らず認識系は結局学習の結果ですから、普段の定位認識の意識度で聞こえ方も変わってきます。「正解」を想像して学習しない限り、脳のネットワークが形成されないのです。
そういう私も開発側なので、デモを聞いてもらうときはいつもドキドキします。
空間音響系の開発側は、「他の人にはけっして同じようには聞こえない」ということを念頭に置く必要があります。
否定的なことばかり書いているかもしれませんが、なにもバーチャルサラウンド系を否定してるわけではありません。むしろそう聞こえたいのです！手軽に3Dオーディオを楽しみたいのです！
でも実際は難しい・・。
応用先
ゲームと音楽では条件がかなり違います。
ゲーム
－オーディオ専用プロセッサの存在
単にCPU/GPUパワーアップだと、まずは映像処理に使われて、「音は余った分でやってね」になりがちです
－対象が、効果を感じやすい雑音系の効果音が多い
－対象はオブジェクト単位
－元と多少音質が変わっても許される
－敵の位置が正確に分かる、没入感が得られるなど、ユーザーメリットが大きい
と、好条件が並ぶのに対し音楽では、
音楽
－2chソースがメイン
－音質が変わるのは好まれない
－そもそも、360度立体音響で音楽を聴きたいという需要がどれだけあるのか？
（もちろんそういうジャンルがあっても良いとは思いますが）
－個人差が大きいため、制作側で処理を掛けられない
と、条件が厳しめです。
Appleの空間オーディオも、少なくとも現段階では、映画・ゲーム・VR/AR等、映像との組み合わせで考えているのではないでしょうか。音楽では逆に、頭動かしても音像が正面から動かないって、イヤホンの意味がね・・。
ヘッドトラッキングして音像も動かすと、HRTFに誤差があっても定位認識率が向上することは知られているので、没入感は得られやすいと思いますが。
MATLABでHRTFをいじろう
前置きが長くなりましたが、実際にHRTFをMATLABで扱ってみましょう。
SOFA(Spatially_Oriented_Format_for_Acoustics)
数年前にAESでHRTF等空間音響用共通フォーマット SOFA が規格化され、各プラットフォーム向け読み込みツールも提供されています。
中身は、サンプリング周波数等の測定条件、各インパルス応答、それに対応するazimuth/elevation等が書かれています（DataType = FIRの場合）。
ここでは、Matlab（/Octave）用の使い方を具体的に見ていきましょう。
－Matlab/Octave API をDLし、任意のフォルダーに解凍します
－MATLABの環境->パスの設定で、API_MO フォルダーをパスに追加します
－SOFAフォーマットのHRTFデータセットを入手
SOFAフォーマットのHRTFデータは例えば以下で公開されています。
・ARI
・RIEC(東北大学先端音情報システム研究室）
R2023bで、SOFA ファイルを直接読み書きできる sofaread
/ sofawrite 関数が追加され、外部 API を使用する必要がなくなりました。
それに伴い、Audio Toolbox に付属する ARI のサンプル HRTF データも、従来の "ReferenceHRTF.mat" に加え、SOFA形式の "ReferenceHRTF.sofa" が追加されています。
ただし外部 API とは返り値の構造が異なるので、sofaread を使う場合は SOFAstart は不要ですが、
SOFAload → sofaread
Data.SamplingRate → SamplingRate
Data.IR → Numerator
に置き換えてください。
また、sofaread() で読み込んだ SOFA オブジェクトから interpolateHRTF() で直接 HRTF 補間もできるようになりました。
APIとデータセットが用意できたら、
SOFAstart;
% sofaFileName = '.sofa'; SOFAフォーマットファイル名を指定
hrtf = SOFAload(sofaFileName);として読み込むと、
hrtf.SourcePosition の1列目がazimuth、2列目がelevationとなります。
azimuthは水平面での正面を0として反時計回りの角度、elevationは仰角になります。
例えばRIECのデータセットで azimuth=45、elevation=0 のHRIRが欲しければ、
SOFAstart;
sofaFileName = 'data/RIEC_hrir_subject_001.sofa';
hrtf = SOFAload(sofaFileName);
azimuth = 45; elevation = 0;
idx = find(hrtf.SourcePosition(:,1)==azimuth & hrtf.SourcePosition(:,2)==elevation);
hrirL = squeeze(hrtf.Data.IR(idx, 1, :));
hrirR = squeeze(hrtf.Data.IR(idx, 2, :));とすれば、hrirL/hrirR にそれぞれ、左耳用／右耳用のインパルス応答（IR=フィルター係数）が入ります。
hrtf.Data.SamplingRateにサンプリング周波数が入るので確認しておきましょう。
RIECのは48kHzなので、48kHzサンプリングの音源を用いる必要があります。
例えばこんな感じで実際にHRTFを掛けた音が出ます。
audioIn = audioread('demo\Five_pan_c2.wav');
audioInMono = mean(audioIn,2);
audioOutput = [conv(hrirL, audioInMono) conv(hrirR, audioInMono)];
% sound(audioInMono, hrtf.Data.SamplingRate);
sound(audioOutput, hrtf.Data.SamplingRate);MATLABなら手間いらず
あるいは AudioToolbox があれば、ARIのHRTFデータが一種類だけ、"ReferenceHRTF.mat" として最初から用意されています。これも48kHzサンプリングです。
また、HRTFデータは当然飛び飛びの音源位置で測定されていますし、測定の都合上、各Elevationに対するAzimuthデータ数は必ずしも同じではありません。
さらに、DBによってその粒度も様々です。
このままでは使いにくいので、MATLABでは間の位置のHRTFを補間する関数 "interpolateHRTF()" も用意されています。
デフォルトのHRTFデータを使う
load 'ReferenceHRTF.mat' hrtfData sourcePosition sampleRate
hrtfData = permute(double(hrtfData),[2,3,1]);
sourcePosition = sourcePosition(:,[1,2]);sofaread を使う場合（R2023b 以降）
hrtf = sofaread('ReferenceHRTF.sofa');
hrtfData = hrtf.Numerator;
sourcePosition = hrtf.SourcePosition(:,[1,2]);
これで、sourcePosition の1列目にはAzimuth、2列目にはElevationが入ります。
この中から欲しい音源位置と同じインデックスのデータをhrtfDataから取り出します。
hrtfData(:,1,:) が左耳用、hrtfData(:,2,:) が右耳用です。
Elevation = 0 の 左耳用 IR を waterfall 表示してみましょう。
figure
dataLen = size(hrtfData,3);
waterfall(1:dataLen, sourcePosition(sourcePosition(:,2)==0,1), squeeze(hrtfData(sourcePosition(:,2)==0,1,:)))
xlabel('Sample')
xlim([1 dataLen])
ylabel('Azimuth (deg)')
ylim([0 360])
zlabel('Amplitude (dB)')
Azimuthに粗密があるのが分かります。
補間してみましょう。
load 'ReferenceHRTF.mat' hrtfData sourcePosition sampleRate
hrtfData = permute(double(hrtfData),[2,3,1]);
sourcePosition = sourcePosition(:,[1,2]);
desiredEl = 0;
Fs = sampleRate;
dAdimuth = 2.5;
dAdimuth_view = 45;
figure(1)
waterfall(1:size(hrtfData,3), sourcePosition(sourcePosition(:,2)==0,1), squeeze(hrtfData(sourcePosition(:,2)==0,1,:)))
xlabel('Sample')
xlim([1 256])
ylabel('Azimuth (deg)')
ylim([0 360])
zlabel('Amplitude')
title('HRIRs Elevation=0 (ARI)')
lx = logspace(1,5,512)*pi/Fs/2;
wf = zeros(360/dAdimuth, 512);
fx = lx/(2*pi)*Fs/1000;
leftIRwf = zeros(360/dAdimuth, 256);
rightIRwf = zeros(360/dAdimuth, 256);
cnt = 1;
for desiredAz=0:dAdimuth:360-1
desiredPosition = [desiredAz desiredEl];
interpolatedIR = interpolateHRTF(hrtfData,sourcePosition,desiredPosition,"Algorithm","VBAP");
leftIR = squeeze(interpolatedIR(:,1,:))';
rightIR = squeeze(interpolatedIR(:,2,:))';
h_L = freqz(leftIR,1,lx);
wf(cnt,:) = 20*log10(abs(h_L));
leftIRwf(cnt,:) = leftIR;
rightIRwf(cnt,:) = rightIR;
cnt = cnt + 1;
end
figure(2)
wy = 0:dAdimuth:360-1;
waterfall(1:256, wy, leftIRwf)
xlabel('Sample')
xlim([1 256])
ylabel('Azimuth (deg)')
ylim([0 360])
zlabel('Amplitude')
title('HRIRs Elevation=0 (ARI) interpolated')
figure(3)
% wy = 0:dAdimuth:360-1;
waterfall(fx, wy, wf)
xlabel('Frequency (kHz)')
yticks(0:dAdimuth_view:360-1);
ylabel('Azimuth (deg)')
ylim([0 360])
zlabel('Gain (dB)')
title('HRTFs Elevation=0 (ARI) interpolated')wy の設定場所が違っていたので figure(2) の後に移動しました。
sofaread を使う場合（R2023b 以降）
interpolatedIR = interpolateHRTF(hrtfData,sourcePosition,desiredPosition,"Algorithm","VBAP");
を以下のように書き換えてください。
interpolatedIR = interpolateHRTF(hrtf,desiredPosition,"Algorithm","VBAP");
interpolateHRTF() の引数にある "VBAP（Vector Base Amplitude Panning）"というのは、空間音響ではよく使われる3次元音響パニング法です。
古いリリースではこのVBAP補間モードにバグがあるようですので、最新にしておきましょう。
ついでに補間後の周波数特性も。
MATLABをお持ちの方はマウスでグリグリ動かして眺めてみてください。
任意の音源位置のHRTFを用いてストリーミング処理するには、例えば以下の様にします。256タップ程度のFIRですから、一つの音源に対して並行して複数のHRTF処理を行い、それを滑らかに切り替えたいとかでなければ、FastConvolver 等を使って周波数軸に変換する必要もないでしょう。
load 'ReferenceHRTF.mat' hrtfData sourcePosition sampleRate
hrtfData = permute(double(hrtfData),[2,3,1]);
sourcePosition = sourcePosition(:,[1,2]);
desiredEl = 0; desiredAz = 90;
HRTF_Position = [desiredAz desiredEl];
interpolatedIR = interpolateHRTF(hrtfData, sourcePosition, HRTF_Position, "Algorithm","VBAP");
leftIR = squeeze(interpolatedIR(:,1,:))';
rightIR = squeeze(interpolatedIR(:,2,:))';
HR_fir_L = dsp.FIRFilter(leftIR);
HR_fir_R = dsp.FIRFilter(rightIR);
% inputFilename = ''; 48kHzサンプリングの音声ファイル名を指定
hAudioSource = dsp.AudioFileReader(inputFilename);
hAudioOut = audioDeviceWriter('SampleRate', 48000);
while ~isDone(hAudioSource)
in = hAudioSource();
m = mean(in,2);
out(:,1) = HR_fir_L(m);
out(:,2) = HR_fir_R(m);
hAudioOut(out);
end
release(hAudioSource);
release(hAudioOut);VSTプラグイン化
最後に、MATLAB標準のHRTFを使った音源パニングを、通常のPowerパニング（LRの二乗和が1）と比べるVSTを作ってみましょう。
以前の記事をご覧になれば、1時間も掛からず作れると思います。
ソースがステレオの場合は内部でモノラルに変換しています。
このHRTFをそのまま掛けると音量がかなり小さくなってしまうので、元とあまり変わらないよう、16倍してあります。
HRIRが48kHzサンプリングデータなので、ソースも48kHzのものを使用してください。
ちょっと試すだけなら、44.1kHzでも構わないと思いますが。
－ソースコード
classdef HRTF_Panning < audioPlugin % Inherit from audioPlugin
properties
% GUI parameters
Angle = 0, Elevation = 0;
BYPASS= false;
PAN = 'HRTF';
end
properties (Constant)
% set GUI properties and layout
PluginInterface = audioPluginInterface( ...
'PluginName','HRTF Rotation',...
'VendorName', 'AudiiSion Sound Lab. LLC.', ...
'VendorVersion', '0.0.1', ...
'UniqueId', 'hiro',...
audioPluginParameter('BYPASS', 'Mapping', {'enum','OFF','ON'}, 'Layout', [1 1]), ...
audioPluginParameter('PAN','DisplayName','PAN', 'Mapping',{'enum','HRTF','PAN'}, 'DisplayNameLocation','none','Style','Dropdown', 'Layout', [2 1]), ...
audioPluginParameter('Angle', 'Mapping',{'int', -180, 180},'Style','rotaryknob', 'Layout', [1 2]), ...
audioPluginParameter('Elevation', 'Mapping',{'int', -45, 90},'Style','rotaryknob', 'Layout', [1 3]), ...
audioPluginGridLayout('RowHeight', [70 15], ...
'ColumnWidth', [80 80 80], ...
'Padding', [20 20 20 20]) ... % [left, bottom, right, top]
);
end
properties(Access = private)
Azimuth = 0;
hrtfData, sourcePosition;
HR_fir_L, HR_fir_R;
end
methods
% Constructor
function plugin = HRTF_Panning
plugin@audioPlugin;
s = coder.load('ReferenceHRTF.mat', 'hrtfData', 'sourcePosition'); % sampleRate=48kHz
ir = s.hrtfData; pos = s.sourcePosition;
plugin.hrtfData = permute(double(ir),[2,3,1]);
plugin.sourcePosition = pos(:,[1,2]);
plugin.HR_fir_L = dsp.FIRFilter([1, zeros(1,255)]);
plugin.HR_fir_R = dsp.FIRFilter(zeros(1,256));
setHRIR(plugin); % Init
end
function set.Angle(plugin,val)
plugin.Angle = val;
if plugin.Angle > 0
plugin.Azimuth = 360 - plugin.Angle; %#ok
else
plugin.Azimuth = -plugin.Angle; %#ok
end
setHRIR(plugin);
end
function set.Elevation(plugin,val)
plugin.Elevation = val;
setHRIR(plugin);
end
%% main
%------ main process -------
function out = process(plugin,in)
m = mean(in,2);
out = [m, m];
if ~plugin.BYPASS
if strcmp(plugin.PAN,'HRTF')
out(:,1) = plugin.HR_fir_L(m);
out(:,2) = plugin.HR_fir_R(m);
else
out(:,1) = cos(((plugin.Angle+90)/2)/180*pi) * m;
out(:,2) = sin(((plugin.Angle+90)/2)/180*pi) * m;
end
end
end
%---------------------------
function reset(plugin)
reset(plugin.HR_fir_L); reset(plugin.HR_fir_L);
setHRIR(plugin);
end
end
methods (Access = private)
function setHRIR(plugin)
HRTF_Position = [plugin.Azimuth plugin.Elevation];
interpolatedIR = interpolateHRTF(plugin.hrtfData, plugin.sourcePosition, HRTF_Position, "Algorithm","VBAP");
leftIR = squeeze(interpolatedIR(:,1,:))';
rightIR = squeeze(interpolatedIR(:,2,:))';
DCgain = 1/16;
plugin.HR_fir_L.Numerator = leftIR / DCgain;
plugin.HR_fir_R.Numerator = rightIR / DCgain;
end
end % methods end
end－VST2プラグイン（32ビット）
サンプル音源
下の動画が、このVSTを使って作った、通常のPowerパニングとHRTFパニングの比較です。6s毎に、通常 → HRTFが切替ります。前半はL90度、後半はR90度です。
いかがでしょう？
私は、通常のパニングより外に広がりますがシャリシャリキンキンした音になってしまい、少なくともこれで音楽を聴く気にはなれません。音像も上がって聞こえます。
自分のHRTFを使えば、元と変わらない音に聞こえるんですかね？
本当はイヤホン特性も補正する必要があるのですが、その影響がどの程度なのかは私も分かりません。
論文等でも音像位置が正しく認識できるかどうかの話ばかりで、音質についてはあまり書かれていないんですよね。
自分のHRTFで音楽を処理して聞いたことがある方がいらっしゃったら、ぜひどうだったか教えてください。
ではまた。