球形正十二面體聲源測量近場HRTF的誤差*

余光正 陳澤為

(華南理工大學(xué)聲學(xué)研究所，廣東廣州510640)

自由場條件下，點(diǎn)聲源產(chǎn)生的聲波經(jīng)頭部、軀干和耳廓等人體生理結(jié)構(gòu)散射后到達(dá)雙耳，構(gòu)成一個線性時不變的聲學(xué)傳輸系統(tǒng)，并可描述為頭相關(guān)傳輸函數(shù)(HRTF).當(dāng)聲源與頭中心的距離不小于1.0m時，HRTF近似與聲源距離無關(guān)，稱為遠(yuǎn)場HRTF;當(dāng)距離小于1.0 m時，HRTF隨距離變化明顯，稱為近場HRTF.HRTF包含聲源的空間定位信息，因而對雙耳聽覺的研究具有重要意義，并在虛擬聽覺重放中有重要的應(yīng)用［1-2］.

實(shí)驗(yàn)測量是準(zhǔn)確獲取HRTF的主要方法.但在近場HRTF測量中，由于聲源靠近受試者，聲源尺度引起的多重散射和指向性誤差將比遠(yuǎn)場測量時更明顯.為減小聲源引起的誤差，一些研究中已開始采用幾種小尺度聲源［3-5］.這些小尺度聲源雖能減少多重散射和改善無指向特性，但在1 kHz以下的低頻特性變差，不能準(zhǔn)確測量低頻的近場 HRTF，而HRTF的低頻特性對聲源距離定位卻十分重要［6］.

對于一般的揚(yáng)聲器系統(tǒng)，改善低頻特性與減少聲源尺度經(jīng)常是矛盾的.筆者采用脈動球形聲源和球形頭部的組合模型評估了聲源多重散射引起的HRTF測量誤差，發(fā)現(xiàn)當(dāng)以1.0dB的幅度誤差為限，聲源距離不小于0.20 m，頻率上限為20 kHz時，球形聲源等效半徑應(yīng)不大于0.05 m［7］.但理想的脈動聲源(球面波聲源)并不存在，實(shí)際的聲源在高頻具有一定指向性.為近似實(shí)現(xiàn)無指向的球面波聲源，筆者設(shè)計制作了半徑為0.035m的球形正十二面體聲源，其工作頻率范圍為0.35～20.00 kHz，在低于8kHz的頻段近似無指向(±3dB)，較已有聲源有較明顯的改善［8］.

在近場HRTF測量中，聲源多重散射和指向性引起的誤差是共存的.因此，文中提出球形正十二面體聲源和球形頭部組合模型，用球諧函數(shù)多極展開和球基函數(shù)旋轉(zhuǎn)變換方法，計算和分析聲源多重散射和指向性對近場HRTF測量的共同影響，并提出了改進(jìn)辦法，最后用自制球形頭部模型和球形正十二面體聲源的測量結(jié)果對理論計算結(jié)果進(jìn)行了驗(yàn)證.

1 HRTF的計算模型及求解

1.1 剛球模型HRTF的定義

文中所采用的坐標(biāo)系統(tǒng)如圖1所示，聲源中心位置用球坐標(biāo)(r0，α0，β0)或空間位置矢量r0表示.其中聲源距離r0通常應(yīng)大于0.2m;仰角α0的取值范圍為0°～180°;方位角β0的取值范圍為－90°～270°.α0=90°表示水平面;水平面上，β0= － 90°和90°分別表示正左方和正右方.在該坐標(biāo)系，將圖1的球形正十二面體聲源替換為點(diǎn)聲源，即可得到剛球模型［6］精確的 HRTF，簡記為H*.考慮到模型的對稱性，左、右耳受聲點(diǎn)(即EL和ER)的HRTF統(tǒng)一定義為

式中:f表示頻率，P*表示點(diǎn)聲源在左、右耳受聲點(diǎn)處的聲壓，P*0表示點(diǎn)聲源在剛球頭部中心的自由場輻射聲壓.剛球模型HRTF的計算方法可參考文獻(xiàn)［6］.

1.2 計算模型及誤差公式

球形頭部(文中將半徑固定為0.0875 m)和球形正十二面體聲源組成的計算模型如圖1所示.

圖1 剛球頭部模型和球形正十二面體聲源的組合模型Fig.1 Combined model consisting of a rigid spherical head model and a spherical dodecahedral sound source

圖1中，E表示空間任意受聲點(diǎn);在以O(shè)和O'為原點(diǎn)的球坐標(biāo)系中，分別用r和r'表示E點(diǎn)的位置矢量，用(r，α，β)和(r'，α'，β')表示E點(diǎn)的球坐標(biāo); 而(x，y，z)和(x'，y'，z') 分別為對應(yīng)的笛卡爾坐標(biāo).對于圖1的計算模型，其HRTF(記為H)可參考式(1)定義為

式中，P表示球形正十二面體聲源在左、右耳受聲點(diǎn)的聲壓，P0為球形正十二面體聲源在球形頭部中心的自由場輻射聲壓.

文中主要分析和討論由球形正十二面體聲源的散射和指向性共同引起的近場HRTF測量誤差，因此這里先給出HRTF測量的幅度譜誤差定義:

式中:SD表示幅度譜誤差，單位為dB.要求解式(3)的SD，須先求解式(2)的組合模型HRTF，即求解式(2)的雙耳聲壓P和自由場聲壓P0.P0的求解可參考文獻(xiàn)［8］，下面簡要推導(dǎo)雙耳聲壓P的求解公式.

1.3 雙耳聲壓P的求解

在圖1中，任意受聲點(diǎn)E的空間位置矢量用r表示(左、右耳的空間位置矢量rL和rR是其特例).受聲點(diǎn)E的復(fù)數(shù)聲壓P(r，r0，f)應(yīng)滿足 Helmholtz方程:

式中:k=2πf/c，表示波數(shù);c表示空氣中的聲速.球形頭部表面SH(設(shè)為剛性聲學(xué)邊界［9］)和球形正十二面體表面SD(局部脈動單元以外的區(qū)域設(shè)為剛性聲學(xué)邊界)的聲壓分別滿足以下邊界條件:

式中:?/?n表示法向?qū)?shù);j為單位復(fù)數(shù);ω=2πf，為角頻率;ρ0為空氣密度;v為球形正十二面體聲源表面的法向速度，且有

式中:i表示局部脈動單元的編號，i=1，2，…，12;v0表示局部脈動單元法向振幅;θ0表示局部脈動單元的半張角;αi為圖2 所示球坐標(biāo)系(ri，αi，βi)中的仰角.在無窮遠(yuǎn)，聲壓P(r，r0，f)滿足 Sommerfeld輻射條件:

式中，?/?r表示對距離求偏導(dǎo)數(shù).求解聲壓P(r，r0，f)也就是求解式(4)在邊界條件(5)－(8)下的解.一般情況下，該聲壓可分解為3部分的和:

式中:P0(r，r0，f)表示聲源在受聲點(diǎn)E的直達(dá)聲壓;PH(r，r0，f)表示球形頭部的散射聲壓;PD(r，r0，f)表示球形正十二面體聲源的散射聲壓.三者可按球諧函數(shù)展開為表達(dá)式(11) －(12)［10］.其中，式(11)在球坐標(biāo)系(r，α，β)中展開，式(10)和(12)在球坐標(biāo)系(r'，α'，β')中展開(見圖 1).

式中，hl是l階第二類球漢克爾函數(shù)，Ylm是復(fù)球諧函數(shù).求解聲壓P(r，r0，f)轉(zhuǎn)而成為求解式(10)－(12)的展開系數(shù)Dlm、Alm、Blm.實(shí)際計算時，先解出系數(shù)Dlm，Alm和Blm則可通過邊界條件式(5)、(6)聯(lián)立方程組求解［7］.

求解式(10)中的展開系數(shù)Dlm時，分別將12個局部脈動單元的輻射聲壓按式(10)展開并線性疊加.當(dāng)?shù)趇個局部圓形脈動單元的中心不經(jīng)過z'軸時，須引入圖 2 所示的球坐標(biāo)(ri，αi，βi)，它對應(yīng)的笛卡爾坐標(biāo)為(xi，yi，zi)，此時，第i個局部圓形脈動單元的輻射聲壓PDi(ri，f)可表示為［11］

式中，展開系數(shù)Dil的求解公式為

aD為球形正十二面體聲源半徑，h'l是l階第二類球漢克爾函數(shù)的導(dǎo)數(shù)，系數(shù)Cl滿足

Γl是l階的勒讓德多項(xiàng)式，且Γ－1=1.

為了把球坐標(biāo)(ri，αi，βi)下的表達(dá)式(13)用球坐標(biāo)(r'，α'，β')重新表達(dá)，需用到以下的球基函數(shù)旋轉(zhuǎn)變換［12］(示意圖見圖2):

式中，α'i和β'i表示第i個局部脈動單元中心在球坐標(biāo)系(r'，α'，β')中的方位角和仰角，變換系數(shù)T可用下面的公式求解［12］:

圖2 球基函數(shù)旋轉(zhuǎn)變換Fig.2 Rotation transform of the spherical basis functions

利用公式(13)－(17)，可將任意第i個局部脈動單元的輻射聲壓表示在球坐標(biāo)系(r'，α'，β')中，然后對12個局部脈動單元的輻射聲壓求和，便可得到式(10)的展開系數(shù)Dlm:

將Dlm代入式(10)－(12)，可結(jié)合邊界條件求出各項(xiàng)展開系數(shù)，從而求解出雙耳聲壓和HRTF.

2 HRTF的誤差分析

2.1 誤差的算例

為簡單驗(yàn)證計算結(jié)果的準(zhǔn)確性，先對幅度譜誤差SD分步進(jìn)行求解，表達(dá)式如下:

式中，SD1和SD2按下式確定:

H1是忽略式(9)的PD(r，r0，f)項(xiàng)后算得的 HRTF，因而第1步求解得到的SD1只反映聲源多重散射的影響，第2步求解的SD2只反映聲源指向性的影響.作為算例，球形正十二面體聲源半徑aD取為0.035m，脈動單元半張角θ0為23.6°，聲源位于r0=0.2 m、α0=90°、β0=90°處，算得右耳受聲點(diǎn)ER的 HRTF幅度譜誤差 SD1、SD2、SD，如圖3 所示.

由圖3可知:SD1在±0.5 dB的范圍內(nèi)，SD2在±2.5dB的范圍;SD恰好是SD1和SD2的疊加.且SD1和SD2的結(jié)果與文獻(xiàn)［7］和［8］的計算結(jié)果基本一致，說明計算結(jié)果是準(zhǔn)確的.

圖3 右耳HRTF的 SD1、SD2和SDFig.3 SD1，SD2and SD of HRTF of right ear

2.2 總體誤差與聲源改進(jìn)

隨著聲源靠近球體，聲源散射誤差將變得明顯［7］，因而這里只取0.2 m(一般作為分析近場HRTF的最小聲源距離［2］)進(jìn)行分析.由圖1所示計算模型的對稱性，只給出α0=90°(水平面)、聲源方位角β0在－90°～90°(前半周)時的右耳 HRTF幅度譜誤差，見圖4.

aD=0.035m和θ0=23.6°時的計算結(jié)果如圖4(a)所示.可見，明顯的譜誤差出現(xiàn)在約8 kHz以上;隨著聲源偏離受聲耳(β0偏離90°)，HRTF的幅度譜誤差略有增大，但出現(xiàn)誤差的頻段基本保持一致.

為了改進(jìn)現(xiàn)有的聲源，另計算分析了兩種情況下的聲源誤差.圖4(b)是減小聲源尺寸至aD=0.025m且保持脈動單元半張角不變時的計算結(jié)果;圖4(c)給出了脈動單元的半張角增大到θ0=30°(接近理想的最大值31.7°)時的計算結(jié)果.由圖中可見，隨著聲源尺寸減小和局部脈動單元半張角增大，球形正十二面體聲源多重散射和指向性引起的誤差幅度減小，且出現(xiàn)誤差的頻率升高，這對近場HRTF的測量是有利的.如圖4(b)所示，當(dāng)聲源半徑為0.025 m、局部脈動單元半張角為23.6°時，較明顯的誤差(一般指大于1 dB)出現(xiàn)在約12 kHz以上，除個別方向( －90°≤β0≤ －50°)和頻率(f≥10kHz)，誤差幅度可控制在 ±3.5 dB.事實(shí)上，在10kHz以上的頻率，HRTF所包含的譜特性對定位的貢獻(xiàn)已逐漸減少，而各種因素(與聲源無關(guān))導(dǎo)致測量HRTF的幅度譜誤差較大(通常在10 dB左右).因此，可認(rèn)為圖4(b)的聲源引起的誤差已基本滿足近場HRTF的測量要求，為進(jìn)一步改進(jìn)聲源提供了理論依據(jù).

圖4 r0=0.2m，α0=90°時右耳HRTF的SDFig.4 SD of HRTF of right ear when r0=0.2m and α0=90°

3 測量驗(yàn)證

為進(jìn)一步驗(yàn)證上述計算結(jié)果的合理性，使用實(shí)際的球形正十二面體聲源(半徑aD約0.035m，脈動單元半張角取θ0=23.6°)和木質(zhì)剛球模型(半徑約0.0875 m)進(jìn)行測量驗(yàn)證，如圖5所示(聲源位于r0=0.2m、α0=90°、β0=90°處).HRTF 的測量環(huán)境和方法與文獻(xiàn)［8］相同.另一方面，按照圖5實(shí)際模型的尺寸，計算圖1所示組合模型條件下的HRTF.測量和計算的HRTF結(jié)果見圖6.由圖6可見，在0.35～10kHz的頻段，測量值與計算值基本一致，聲源指向性和散射的總誤差約在±2.5 dB以內(nèi).測量值在0.35Hz以下明顯衰減，這是由聲源的頻率響應(yīng)特性所決定的.其它譜細(xì)節(jié)差異是由實(shí)驗(yàn)誤差(主要是實(shí)際制作的聲源和球體與理想模型的差異)引起的.上述實(shí)驗(yàn)結(jié)果進(jìn)一步驗(yàn)證了計算結(jié)果的合理性.

圖5 剛球模型以及球形正十二面體聲源Fig.5 Rigid spherical model and spherical dodecahedral sound source

圖6 測量和計算的右耳HRTFFig.6 Measured and calculated HRTF of right ear

4 結(jié)語

文中提出了球形正十二面體聲源和球形頭部的組合計算模型，利用球諧函數(shù)多極展開和球基函數(shù)旋轉(zhuǎn)變換的方法，計算和分析了近場HRTF測量中由球形正十二面體聲源指向性和多重散射共同引起的HRTF幅度譜誤差.計算結(jié)果的準(zhǔn)確性通過兩步求解的辦法得到了驗(yàn)證，計算結(jié)果的合理性用實(shí)際制作的球形正十二面體聲源和木質(zhì)球體進(jìn)行了實(shí)驗(yàn)驗(yàn)證.文中還在不同方向全面分析了聲源對近場HRTF測量的影響，結(jié)果表明，隨著聲源尺寸的減小和局部脈動單元半張角的增大，球形正十二面體聲源多重散射和指向性引起的誤差幅度減小，且有向高頻移動的趨勢.當(dāng)聲源距離不小于0.2 m、聲源半徑不大于0.025 m、局部脈動單元半張角不小于23.6°時，較明顯的誤差出現(xiàn)在約12 kHz以上，且誤差幅度基本在±3.5 dB，可滿足近場HRTF的測量要求.文中結(jié)果為進(jìn)一步改進(jìn)適合于近場HRTF測量的球形正十二面體聲源提供了理論依據(jù).

［1］謝菠蓀.頭相關(guān)傳輸函數(shù)與虛擬聽覺［M］.北京:國防工業(yè)出版社，2008:1-33.

［2］余光正，謝菠蓀.近場頭相關(guān)傳輸函數(shù)及其應(yīng)用［J］.電聲技術(shù)，2007，31(7):45-50.Yu Guang-zheng，Xie Bo-sun.Head-related transfer function for nearby sources and its application［J］.Audio Engineering，2007，31(7):45-50.

［3］Nishino T，Hosoe S，Takeda K，et al.Measurement of the head related transfer function using the spark noise［C］∥The 18th International Congress on Acoustics.Kyoto:Science Council of Japan，2004:1437-1438.

［4］Hosoe S，Nishino T，Itou K，et al.Development of micrododecahedral for measuring head-related transfer functions in the proximal region［C］∥Proceedings of the 2006 IEEE International Conference on Acoustics，Speech，and Signal Proceedings.Toulouse:Centre de Congrès Pierre Baudis，2006:329-332.

［5］Brungart D S，Rabinowitz W M.Auditory localization of nearby sources.head-related transfer functions［J］.Journal of the Acoustic Society of America，1999，106(3):1465-1479.

［6］Duda R O，Martens W L.Range dependence of the response of a spherical head model［J］.Journal of the Acoustic Society of America，1998，104(5):3048-3058.

［7］Yu G Z，Xie B S，Rao D.Effect of sound source scattering on measurement of near-field head-related transfer function ［J］.Chinese Physics Letters，2008，25(8):2926-2929.

［8］Yu G Z，Xie B S，Rao D.Directivity analysis on spherical regular polyhedron sound source used for near-field HRTF measurements［J］.Chinese Physics Letters，2010，27(12):124302.

［9］余光正，謝菠蓀，饒丹.軀干散射對近場HRTF影響的微擾法分析［J］.華南理工大學(xué)學(xué)報:自然科學(xué)版，2010，38(3):143-147.Yu Guang-zheng，Xie Bo-sun，Rao Dan.Analysis of effect of torso scattering on near-field head-related transfer functions by using perturbation method ［J］.Journal of South China University of Technology:Natural Science Edition，2010，38(3):143-147.

［10］Gumerov N A，Duraiswami R.Computation of scattering from N spheres using multipole reexpansion［J］.Journal of the Acoustic Society of America，2002，112(6):2688-2701.

［11］張海瀾.理論聲學(xué)［M］.北京:高等教育出版社，2007:244-248.

［12］Gumerov N A，Duraiswami R.Fast multipole methods for the Helmholtz equation in three dimensions［M］.Oxford:Elsevier，2004:122-124.