見出し画像

MATLABでHRTF~頭部伝達関数とは?~聴覚の仕組みは解明されていない~

前書き

さて最近、今まで一般の場ではほとんど聞かれることのなかった、「HRTF(頭部伝達関数)」という単語がちらほらツイッター等でも見かけるようになりました。
そう、あのゲーム機の影響です。

私はHRTFが専門というわけではなく、自分のHRTFを測定したこともありませんが、何度か試してみたことはあるので、簡単な使い方を解説してみようかと思います。

HRTFとは

HRTF自体は1970年代前半からある古い技術です。

音は音源の位置によって、左右の耳に届く時間や音量がわずかに違ったり、胴体で反射したり、耳介のくぼみ等で共鳴したり小さく反射したり、横方向からであれば、逆側の耳へは頭で遮られたり、頭を回り込んだり、様々な経路を通って、その人の体の影響を受けてから耳に届きます。
人間は(他の動物も)、これらの影響の違いを学習して、音源方向を推定していると言われています。
「これらの影響」を音源から耳への経路の特性として表したのが、頭部伝達関数(Head Related Transfer Function)です。

この特性はかなり複雑で、体の影響を受けているので個人個人全く違っているという特徴があります。

ですので、もし他の人のHRTFが掛かった音を聞くと、音源方向を正しく認識できなかったり、音質が大きく変わってしまったりします

また、耳介を粘土等で形状を変えると、最初は音源方向が曖昧になりますが一ヶ月もすれば正しい定位認識ができるようになり、元に戻しても定位認識はできる、と言う報告もあります。その論文では、「上書きではなく、新しい言語を覚えるように両方のHRTFを認識できるようになる」とあります。

HRTFの周波数特性の違いの例を、以下にグラフで示します。

画像1

音源位置(水平角度)による違い

画像2

人による違い

画像3

同一人物の左耳、右耳による違い


実はステレオスピーカーで音を聞く場合は、「ちゃんと」セッティングしてあげれば、普通のステレオ音源でもかなり立体的な音像として聴くことができます。
これは、「自分の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(東北大学先端音情報システム研究室)

APIとデータセットが用意できたら、

SOFAstart;
% sofaFileName = '.sofa';  SOFAフォーマットファイル名を指定
hrtf = SOFAload(sofaFileName);

として読み込むと、
hrtf.SourcePosition の1列目がazimuth、2列目がelevationとなります。
azimuthは水平面での正面を0として反時計回りの角度、elevationは仰角になります。

画像9

画像9


例えば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]);

これで、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)')

画像8

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)
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')

画像8

interpolateHRTF() の引数にある "VBAP(Vector Base Amplitude Panning)"というのは、空間音響ではよく使われる3次元音響パニング法です。
古いリリースではこのVBAP補間モードにバグがあるようですので、最新にしておきましょう。

ついでに補間後の周波数特性も。
MATLABをお持ちの方はマウスでグリグリ動かして眺めてみてください。

画像8


任意の音源位置の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度です。

画像9

VSTパラメータオートメーション設定DAW画面


いかがでしょう?

私は、通常のパニングより外に広がりますがシャリシャリキンキンした音になってしまい、少なくともこれで音楽を聴く気にはなれません。音像も上がって聞こえます。

自分のHRTFを使えば、元と変わらない音に聞こえるんですかね?
本当はイヤホン特性も補正する必要があるのですが、その影響がどの程度なのかは私も分かりません。
論文等でも音像位置が正しく認識できるかどうかの話ばかりで、音質についてはあまり書かれていないんですよね。

自分のHRTFで音楽を処理して聞いたことがある方がいらっしゃったら、ぜひどうだったか教えてください。

ではまた。





この記事が気に入ったら、サポートをしてみませんか?
気軽にクリエイターの支援と、記事のオススメができます!
3
メーカーの研究所でオーディオ/ビデオのデジタル信号処理アルゴリズム関連 29年 → オーディオ系ベンチャー 5年 → AudiiSion Sound Lab. LLC. 創業(2020/10~) MATLAB暦は多分15年以上? 趣味は楽器演奏・写真