尊重原创，请勿转载! 作者：图林根の烤肠，如有纰漏欢迎指出。

日期：2015年12月31日

文章来源：www.mahong.me/archives/123

2020年6月注：因为年代是在久远，手头目前在忙于Unity和AR项目的事情，待时间赋予会抽空写一篇通过Python实现的内容，还请见谅～

讲完了上一篇，现在可以开始实际演练。具体代码都可以从我的GitHub上面的HRTF文件夹下下载。

如果无法下载，可以移步百度网盘：

链接: https://pan.baidu.com/s/1dslhSYPa3bqSkioPjJtwpw 密码:lyi6

首先主要介绍一下这些文件的功能：

doc 文件夹里面是关于CIPIC HRTF的文档说明，具体相关的资料也可以参考官方网站，个人建议下载官方的那个全部Matlab文档（MATLAB version），总计170 Mb左右, 在下载完成后可以运行里面的showdata.m 文件，从而查看各个位置相关的频谱信息。
audio_with_brir.m 该函数是主要算法，描述了音频文件如何与 HRTF数据进行卷积，从而实现3D音效
hrir_final.mat 是CIPCI KEMAR 21号（Subject_021）的HRTF数据，该数据采用了大耳廓进行采集。
load_CIPIC_HRIR.m 的主要功能是读取 hrir_final.mat 的左右耳HRIR数据，然后在用于audio_with_brir.m文件进行信号卷积
nokia.wav 这个不用多说了，一段测试音频。诺基亚的经典铃声
main.m 为代码的运行的主文件

下面开始分别对各个文件进行解释说明：

a) 读取HRTF函数

从CIPIC的官方文档来看，该HRTF数据的空间坐标系为：

Azimut（方位角，及水平范围内）取值范围【-80 -65 -55 -45 ：5：45 55 65 80】一共25个点，其中从-45°到45°梯度为5°。以人眼为参照物来看，-90°为左耳方向，+90°为右耳方向。

Elevation（立体角？垂直范围内）的取值范围为【 -45 + 5.625 * （0：49）】一共50个采样点。说得有些抽象，看下面一张图便一路了然。

图a从左往右数有25个点（25列，可以想像把手臂伸直，指向你的正前方，然后左右移动手臂画出平行于眼前的一个弧形，这个扇形平面由25个点组成），与之对应的即为Azimut，图b对应的为Elevation，一共50个采样点(同理手臂伸直，上下摆动组成另一个扇形平面，这个平面一共50横行，从B图中右面逆时针数即可)。

接下我们来看HRTF的数据文件，双击通过Matlab调用hrir_final.mat文件后会发现，hrir_l 与 hrir_r 就是与之对应的左右耳HRIR数据。但是随之来的问题就是hrir_l的文件size是 25 x 50 x 200 ，这是因为HRTF的数据是一个三维空间数据，其中的25表示的时Azimut（水平面 Figure1-a），50为Elevation（垂直面Figure1-b），200为该点的数据。假如我们需要取一个平行于头部正前方和正后方的一个圆形区域，由图b的紫色横线可知，脸部正前方的Elevation 对应的点的数值顺序为9（Figure1-b，逆时针数，紫色的线为第9），请不要忘记Azimut的取值范围为左耳-80°到右耳+80°的一个扇形区域，所以这时候我们得到的数据只是头部前面的半扇形区域。（为了更好地理解，可以把双手分别向左右耳的方向伸展为80°，这时候双手围成的一个水平扇形区域就是图a所示的区域，这时候再把双手一同上下移动，就会构成如图b中一样的一个三维空间）这时候我们把紫线延伸到图a可以看到，头部正后方Elevation 为41的区域，正是我们想要的后方区域。明白了上面的原理，再来看代码第22行

out = [squeeze(hrir_l(:,front,:));flip(squeeze(hrir_l(:,back,:)),1)];      %（load_CIPIC_HRIR.m）

该行代码的意义为：先取头部正前方Elevation 为 front 的25个水平面上的点，再取脸部后方Elevation为back的25个水平面上面的点翻转后在合并，就组成了一个由50个点组成的水平面。至于说为什么要翻转后平面，因为hrtf的取值是按照顺时针采样的，而后半部数据是前半部数据的镜像（不知道这么说算解释么 -_- !），不翻转数据直接进行卷积后会发现声音是从-80->0 ->80 ->-80->-180->80的顺序播放。

b) 卷积前的操作

得到了HRIR数据后我们就可以开始正式进行信号的卷积了，这里我假设您已经阅读完我前一篇的文章并以理解其原理。

首先来看 audio_with_brir.m 函数的33到38行，

if remainder > 0 % if the remainder not equal to zero,then calculate the padding zeros
   padding_zero = points_number - remainder;

   % padding zeros with audio file
   audioIn = [audioIn zeros(1,padding_zero) ];
end

这段代码的意义是为了进行给音频末尾补零，以便能够让输入音频平均分成N等分。算法就是首按照输入的BRIR数据有多少个纵列值（BRIR必须包括左右耳的数据，并且数值为纵列）来进行音频的N等分（这里为points_number）。然后查看音频的总长度除以需要N等分的余数来进行补零。还是举个简单的例子，11/4 余数为 3，则需要补零的个数为4-3 = 1

接下来就是设计窗函数来进行淡入淡出，以及高通与低通的滤波器，这里的截止频率分别为20kHz与20Hz。

第55到64行，L和R分别为左右耳

for i = 1 : points_number % because brir has only 6 points in Azimut   {300, 330, 0, 60, 120, 180}

    % faltung    
    L = filter(brir(:,2*i-1)',1,audioIn(1,step*(i-1)+1:step*i));
    R = filter(brir(:,2*i)',1,audioIn(1,step*(i-1)+1:step*i));

    segment_L(i,:) = L;
    segment_R(i,:) = R;

end

即为音频分割后的segment与BRIR数据进行卷积。

然后是67到78行：为了代码方便，这里的把音频的片段1和片段最后一段也进行了淡入淡出效果，虽然我们不需要第一段开头的淡入和最后一段的淡出，但是加上后也无妨。

for n = 1 : points_number

    % the last segment need only fade in, so we make that later, first make
    % fade in and fade out except the last segment
    % add fade in for each segment at the beginning
    segment_L(n,1:window_N/2) = segment_L(n,1:window_N/2) .* window(1:window_N/2);  
    segment_R(n,1:window_N/2) = segment_R(n,1:window_N/2) .* window(1:window_N/2); 

    %add fade out for each segment at the end
    segment_L(n,(end - window_N/2 + 1):end) = segment_L(n,(end - window_N/2 + 1):end) .* window((window_N/2+1):end); 
    segment_R(n,(end - window_N/2 + 1):end) = segment_R(n,(end - window_N/2 + 1):end) .* window((window_N/2+1):end);    
end

即为信号的淡入淡出操作。

86和87行为第一个segment的操作：

left = [segment_L(1,1:end - (window_N/2)) segment_L(1,(end - window_N/2 +1):end) + segment_L(2,1:(window_N/2))];  
right = [segment_R(1,1:end - (window_N/2)) segment_R(1,(end - window_N/2 + 1):end) + segment_R(2,1:(window_N/2))];

93到96为直到倒数第二个segment的合并

for n = 2 : points_number - 1 % except the last segment
    left = [left segment_L(n,(window_N/2+1):end-(window_N/2)) segment_L(n,(end - window_N/2 + 1):end) + segment_L(n+1,1:(window_N/2))];
    right = [right segment_R(n,(window_N/2+1):end-(window_N/2)) segment_R(n,(end - window_N/2 + 1):end) + segment_R(n+1,1:(window_N/2))];
end

100到120行为最后一段segment的合并

n = points_number;
left = [left segment_L(n,(window_N/2 +1):end)];
right = [right segment_R(n,(window_N/2 + 1):end)];

余下的代码就是音频的归一化与滤波。

c) 信号卷积

明白了上面两个函数的功能后我们可以开始正式进行信号的卷积操作，这时候来看main文件，

15,16行调用了load_CIPIC_HRIR函数，分别读取了Elevation为9和41的数值

hrir_l = load_CIPIC_HRIR(hrir_fn,front,back,'left');
hrir_r = load_CIPIC_HRIR(hrir_fn,front,back,'right');

22到26行把所得到的所有HRIR数据分别按照左耳1，右耳1，左耳2，右耳2等等整合到了一起。

hrir = [];
for i = 1:size(hrir_l,2)
    hrir(:,i*2-1) = hrir_l(:,i);
    hrir(:,i*2) =  hrir_r(:,i);
end

最后就是audio_with_brir函数，输入参数为音频文件与最后的HRIR。

然后呢？没有然后了，赶紧来听听诺基亚的环绕音效吧。

以上就是最近几天来关于毕设HRTF的一些编程总结，如果想采用其它环绕或者有个性的方法，只需要读取不同空间位置的HRIR数据即可。如仍有其它问题的话，欢迎邮件咨询，

最后祝愿各为2016新年快乐~

图林根の烤肠

图林根の烤肠

小时候喝三鹿，长大用百度

HRTF 3D 音效 Matlab实现

a) 读取HRTF函数

b) 卷积前的操作

c) 信号卷积