CN102664017A - 一种3d音频质量客观评价方法 - Google Patents

Abstract

本发明是一种3D音频质量客观评价方法，经过双耳合成后，分别提取多声道信号的音色感知特征参数、水平声场感知特征参数和高度声场感知特征参数这三种底层声学特征参数，最后通过3D声场心理认知模型的处理以达到拟合主观评价的结果，从而实现客观评价方法与主观评价方法的最大相关。本方法考虑了双耳对三维空间内的声源信息的音色及声场的感知特征参数，客观评价结果更接近三维感知声场的主观评价得分。

Description

一种3D音频质量客观评价方法

技术领域

本发明涉及语音音频质量客观评测领域，尤其涉及3D音频和空间音频质量客观评价方法优化。

背景技术

人类的听觉中，除了对声音的响度、音调和音色等主观属性的感觉外，还包括对声音的空间听觉，也就是对声音空间属性或特征的主观感觉，3D音频就是一种强调了人耳对声音的空间属性主观感受的音频信号。随着3D音频在多媒体领域和娱乐产业中的广泛应用，如何评价3D音频设备的听音感受和音频质量就成了一个亟待关注的研究点。

目前已经有许多音频系统性能评价方法，大致可分为主观和客观评价两大类。ITU组织针对二维环绕声效评价早在1996年和2003年发布了ITU-RBS.1116-1和BS.1534-1两种宽带音频质量主观评价标准。随后，NHK针对其22.2多声道系统在2006年提出了基于语义区分的三维感知声场主观评价方法。

由于主观评价是以人作为主体直接参与评价，其评价结果相对更直接、更合理。传统的音频质量主观评价方法比较典型的有平均意见得分(Mean Opinion Score：MOS)评价方法，MOS评价方法是用五个等级来评价音频质量等级：优(5分)、良(4分)、一般(3分)、差(2分)、坏(1分)。但是主观评价的主要问题是耗时费力，代价大。这使得研究者寻找一种通过建立预测人感知音频质量的模型和方法进行音频质量的客观评价，其结果应该与主观听觉测试越接近越好。

目前国际主流上的音频质量客观评价方法是ITU-R BS.1387-1(PEAQ)，PEAQ客观评价方法是ITU综合已有的六种评价方法：Disturbance Index(DIX)，Noise-to-Mask Ratio(NMR)，Perceptual Audio Quality Measure(PAQM)，Perceptual Evaluation(PERCEVAL)，PerceptualObjective Measure(POM)and The Toolbox Approach后提出的一项针对单声道的音频信号的客观评价方法。根据ITU-RBS.562标准将音频质量等级以5级递减表示：无法感知(5分)、可感知但不恼人(4分)、轻微恼人的(3分)、恼人的(2分)、非常恼人的(1分)。则通过参考信号和待测信号评分相减得到客观评价得分ODG，分值范围为0到-4之间。PEAQ在评价音频质量的音色感知声学属性上与主观评价相近，但是PEAQ只能对单声道信号进行客观评价，不能表征3D音频的空间感知特征参数。

发明内容

本发明的目的是提供3D音频质量客观评价技术方案，使得客观评价的结果更符合人对3D音频的主观感知效果。

基于以上目的，本发明提供一种3D音频质量客观评价方法，包含以下步骤：

步骤1，将多声道信号的参考序列下混得到参考信号的双耳合成信号，将多声道信号的测试序列下混得到待测信号的双耳合成信号；

步骤2，根据步骤1所得参考信号的双耳合成信号和待测信号的双耳合成信号，计算得到待测信号的音色感知特征参数、水平声场感知特征参数、高度感知特征参数；

步骤3，随机选取若干组步骤2所得待测信号的音色感知特征参数、水平声场感知特征参数、高度感知特征参数作为3D声场心理认知模型的输入，以预设的相应主观评价分数作为3D声场心理认知模型的输出，经过训练得到3D声场心理认知模型，将步骤2所得其他待测信号的音色感知特征参数、水平声场感知特征参数、高度感知特征参数输入训练后的3D声场心理认知模型，得到音频质量客观评价结果。

而且，在步骤1中，采用头相关传输函数实现下混，得到双耳合成信号。

而且，采用ITU-R BS.1387-1标准的感知音频质量客观评价模块实现计算待测信号的音色感知特征参数，包括如下步骤，

步骤1-1，通过FFT变换将双耳合成信号变换为频域信号，对频域信号计算平方幅值得到能量信号，模拟人外耳和中耳听觉效应对能量信号进行滤波处理，得到滤波后的能量信号；

步骤1-2，将滤波后的能量信号进行临界带分组，得到临界带分组信号；所述临界带分组处理模拟了人耳基底膜对不同频带信号的响应差异；

步骤1-3，根据心理声学理论对临界带分组信号进行临界带扩展得到激励样本，并计算得到调制样本；

步骤1-4，对激励样本经过时域和频域平滑处理，然后通过电平和模式调整，得到频谱适应样本，并计算掩蔽阈值；

步骤1-5，根据上述四个步骤的计算结果计算每帧音色感知特征参数值，最后通过平均计算得到待测信号的音色感知特征参数。

而且，所述待测信号的水平声场感知特征参数包括双耳声级差异、双耳时间差异和双耳互相关系数，计算水平声场感知特征参数包括如下步骤，

步骤2-1，通过FFT变换将双耳合成信号变换为频域信号，将参考信号和待测信号的频域信号通过一个Gammatone滤波器组划分为24个频段的巴克带，模拟外耳和中耳滤波对划分的巴克带进行加权，得到参考信号的激励模式和待测信号的激励模式；

步骤2-2，计算参考信号的激励模式和待测信号的激励模式的归一化互相关函数，得到参考信号和待测信号之间的双耳时间差异、参考信号和待测信号之间的双耳互相关系数，通过计算参考信号的激励模式和待测信号的激励模式的声强比值的对数值，得到参考信号和待测信号之间的双耳声级差异。

而且，所述待测信号的高度感知特征参数包括高度角，计算高度感知特征参数包括如下步骤，

步骤3-1，分离待测信号的双耳合成信号的左右声道，并分别进行快速傅里叶变换，得到待测信号的频谱信号；

步骤3-2，将待测信号的频谱信号通过一个高斯滤波器滤波，得到滤波后的频谱信号；

步骤3-3，在待测信号的滤波后的频谱信号中检测频谱的波峰和波谷，将出现频谱的波峰波谷的位置和头相关传输函数数据库中的不同高度角的波峰波谷位置进行对比，估计得到待测信号的高度角。

而且，所述3D声场心理认知模型使用BP神经网络，在3D声场心理认知模型中BP神经网络使用最速下降法。

本发明技术方案评价了3D音频质量多个属性，包括整体音频质量(BAQ)，音色保真度(TF)，3维前方空间保真度(3DFSF)，3维环绕空间保真度(3DSSF)和3维高度空间保真度(3DESF)，其中通过音色感知模型的PEAQ测度的MOVs值用来表征整体音频质量(GAQ)和音色保真度(TF)，而通过水平声场感知模型计算的水平声场感知特征参数双耳互相关系数(Inter-Aural Cross-Correlation：IACC)，双耳时间差异(Interaural Time Difference：ITD)和双耳声级差异(Interaural Time Difference：ITD)作为表征3维前方空间保真度(3DFSF)和3维环绕空间保真度(3DSSF)的评测参数，通过高度感知模型计算的高度角θ作为评价3维高度空间保真度(3DESF)的评价参数。

附图说明

图1是本发明实施例的流程图。

图2是本发明实施例的音色感知模型框图。

图3是本发明实施例的水平声场感知模型结构框图。

图4是本发明实施例的高度感知模型结构框图。

图5是本发明实施例的3D声场心理认知模型结构框图。

具体实施方式

以下结合附图和实施例详细说明本发明技术方案。

本发明实施例提供的方法包括以下步骤，流程图参见图1：

步骤1，将多声道信号的参考序列下混得到参考信号的双耳合成信号，将多声道信号的测试序列下混得到待测信号的双耳合成信号。可以采用头相关传输函数实现下混，得到双耳合成信号。

实施例在步骤1将三维重建的多声道信号的参考序列和测试序列经过双耳信号模拟发生器，双耳信号表示各声道的输入与对应扬声器位置的头相关冲击响应(Head Related ImpulseResponses：HRIRs)卷积模拟生成对应双耳接收的两路信号。设N为声道总数，参考和待测信号双耳信号合成如式(1)所示。

$\left(\begin{array}{cc}{\hat{L}}_{\mathrm{TestH}}& {\hat{L}}_{\mathrm{RefH}}\\ {\hat{R}}_{\mathrm{TestH}}& {\hat{R}}_{\mathrm{RefH}}\end{array}\right)=\left(\begin{array}{ccc}{H}_{\mathrm{ch}1\_l}& H_{\mathrm{ch}2\_l}L& H_{\mathrm{chN}\_l}\\ H_{\mathrm{ch}1\_r}& H_{\mathrm{ch}2\_r}& H_{\mathrm{chN}\_r}\end{array}\right)\left(\begin{array}{cc}{\mathrm{ch}1}_{\mathrm{test}}& {\mathrm{ch}1}_{\mathrm{Ref}}\\ {\mathrm{ch}2}_{\mathrm{test}}& {\mathrm{ch}2}_{\mathrm{Ref}}\\ M& M\\ {\mathrm{chN}}_{\mathrm{test}}& {\mathrm{chN}}_{\mathrm{Ref}}\end{array}\right)...\left(1\right)$

其中H_{ch1_l}，

表示第1至N个声道信号抵达左耳的头相关传输函数(HeadRelated Transfer Function：HRTF)，H_{ch1_r}，表示第1至N个声道信号抵达右耳的头相关传输函数。ch1_test，ch2_test...chN_test表示第1至N个声道待测信号的双耳合成信号，ch1_Ref，ch2_Ref...chN_Ref表示第1至N个声道参考信号的双耳合成信号，和

分别表示左右耳的输入信号，上式中头相关冲击响应表示时域而头相关传输函数表示频域。其中，

和分别表示左耳的输入待测信号和输入参考信号、

和

分别表示右耳的输入待测信号输入的参考信号。

不失一般性，在下文中，将上述双耳合成信号ch1_test～chN_test和ch1_ref～chN_ref分别用待测信号的双耳合成信号s_test和参考信号的双耳合成信号s_ref表示。

步骤2，根据步骤1所得参考信号的双耳合成信号和待测信号的双耳合成信号，计算得到待测信号的音色感知特征参数、水平声场感知特征参数、高度感知特征参数。

实施例将计算音色感知特征参数、水平声场感知特征参数、高度感知特征参数的具体过程分别实现为音色感知模型、水平声场感知模型和高度感知模型。用步骤1所得用待测信号的双耳合成信号s_test和参考信号的双耳合成信号s_ref分别作为音色感知模型、水平声场感知模型和高度感知模型的输入，即可输出得到相应的音色感知特征参数、水平声场感知特征参数、高度感知特征参数。

图2说明本发明所述方法音色感知模型框图，参考信号和待测信号的双耳合成信号s_test和s_ref输入音色感知模型，实施例中音色感知模型采用ITU-RBS.1387-1标准的PEAQ模块来评价，输出10个表征音频音色的特征参数MOVs，该模型是基于FFT的人耳模型，具体实现为现有技术。

为便于实施参考起见，提供实施例计算音色感知特征参数的具体步骤如下所述：

步骤1-1，时频变换，即通过FFT变换将双耳合成信号变换为频域信号；然后进行外耳和中耳滤波。

参考信号的双耳合成信号s_ref为时域信号，将参考信号的双耳合成信号s_ref通过短时傅里叶变换，得到参考信号的频域信号S_ref；同样的，待测信号的双耳合成信号s_tesf为时域信号，将待测信号的双耳合成信号s_test通过短时傅里叶变换，得到待测信号的频域信号S_test。实施例中输入的参考信号的时域信号和待测信号的时域信号采样率为F_s＝48kHz，通过短时傅里叶变换得到的频域信号，帧长为N_F＝2048个点。下文为便于用公式描述，将参考信号和待测信号的频域信号都以序列的方式表示为F[k_f，n₁]，其中n₁表示信号帧数，k_f表示信号频率(下同)。

实施例进行的外耳和中耳滤波，包括对频域信号计算平方幅值得到能量信号，模拟人外耳和中耳听觉效应对能量信号进行滤波处理。

外耳和中耳的频率响应通过传输函数来模拟，得到S_ref和S_test的外耳加权FFT输出F_e[k_f，n₁]，即滤波处理后的能量信号。外耳传输函数W(f)/dB如下：

$W\left(f\right)/\mathrm{dB}=-0.6\·3.64\·{\left(\frac{f}{1000}\right)}^{-0.8}+6.5\·e^{0.6\·{(\frac{f}{1000}-3.3)}^2}-{10}^{-3}\·{\left(\frac{f}{1000}\right)}^{3.6}---\left(2\right)$

式中，频率f的单位是HZ，外耳和中耳的传输函数表示了外耳和中耳对声音的衰减和滤波过程。便于计算缘故，用W[k_f]表示窗函数对应于输入信号的序列表示方式，则需要对W(f)进行如下式的转换：

$W\left[k_f\right]=W\left(\frac{{\mathrm{kF}}_s}{N_F}\right)...\left(3\right)$

上式中，

k为采样点下标，则有参考信号和待测信号的外耳加权输出F_e[k_f，n₁]为

$F_e[k_f,n_1]=\left|F\right[k_f,n_1\left]\right|\×{10}^{\frac{W\left[k_f\right\}}{20}}...\left(4\right)$

步骤1-2，临界带分组，包括将步骤1-1所得滤波后的能量信号进行临界带分组，得到临界带分组信号；所述临界带分组处理模拟了人耳基底膜对不同频带信号的响应差异。

本步骤为模拟内耳基底膜的频率特性，频谱被分配到数个互不交叠的子带中，首先将参考信号和待测信号的频域信号建立从频域到Bark域的映射：

$Z_1/\mathrm{Bark}=7\·\mathrm{ar}\sinh \left(\frac{f/\mathrm{Hz}}{650}\right)...\left(5\right)$

其中Z₁为Bark带个数，f为人类听觉可闻范围：80Hz--18000Hz，ar sinh(.)是反正弦函数。Zwicker and Feldtkeller在1967年根据认为听觉特性提出可以将频率划分为活干个临界频带，以及Bark域测度，规定了频率映射到Bark域的规则。所谓临界带分组其实是把DFT能量向量分组对应到临界带，其中能量向量用DFT频域值的平方幅值来表示，这是一种映射关系。

根据公式P_e[k₁，n₁]＝|F_e[k_f，n₁]|²求得外耳加权FFT输出的能量值，这个过程得到每个频带的子带频谱能量值P_e[k₁，n₁]。

实施例还添加内部噪声，给每个子带添加一个频率偏移量PThres，第k₁个子带的频率偏移量记为P_Thres[k₁]，得到音高样本P_p[k₁，n₁]，作为临界带分组信号。计算公式如下：

$P_{\mathrm{Thres}}\left[k_1\right]={10}^{0.4-0.364\×{\left(\frac{f_c\left[k_1\right]}{1000}\right)}^{-0.8}}...\left(6\right)$

$P_p[k_1,n_1]=P_e[k_1,n_1]+P_{\mathrm{Thres}}\left[k_1\right]---\left(7\right)$

其中，f_c[k₁]表示第k₁个子带的中心频率。

步骤1-3，时域和频域扩展，根据心理声学理论对步骤1-2所得临界带分组信号进行临界带扩展得到激励样本，并计算得到调制样本。

频域扩展在Bark域上进行。将音高样本P_p[k₁，n₁]经由一个扩展函数分配至整个听域空间，那么第k₁个子带的能量就是各个子带能力在该子带贡献的加和；引入归一化因子，消除扩展函数在各个子带中引入的增益，得到未拖尾的激励样本E_unsmeared[k₁，n₁]，根据未拖尾的激励样本E_unsmeared[k₁，n₁]计算调制能量和包络大小得到调制样本。为了模拟向前掩蔽，将各个子带的能量在时间上低通滤波器平滑，得到激励样本E[k₁，n₁]。

步骤1-4，激励样本处理，包括对激励样本经过时域和频域平滑处理，然后通过电平和模式调整，得到频谱适应样本，并计算掩蔽阈值。

本步骤要对参考信号和待测信号的激励样本进行电平适应和模式适应处理。为补偿参考信号和待测信号的响度差异和线性失真，待测信号和参考信号的平均电平要适应调整。适应分为电平适应和模式适应两个步骤，最终获得参考信号和待测信号的频谱适应样本。

其次是调制计算部分，首先根据未拖尾的激励样本E_unsmeared[k₁，n₁]计算平均响度

和计算简化响度的差值

计算公式如下：

${\stackrel{\‾}{E}}_{\mathrm{der}}[k_1,n_1]=a\·{\stackrel{\‾}{E}}_{\mathrm{der}}[k_1,n_1-1]+(1-a)\·\frac{F_s}{N_F/2}\·|E{[k_1,n_1]}^{0.3}-E{[k_1,n_1-1]}^{0.3}|...\left(8\right)$

$\stackrel{\‾}{E}[k_1,n_1]=a\·\stackrel{\‾}{E}[k_1,n_1-1]+(1-a)\·E{[k_1,n_1]}^{0.3}...\left(9\right)$

式中

τ表示时间常数，E[k₁，n₁]表示激励样本模式。

然后可得到各子带包络，子带包络表示每一个滤波通道的包络大小，计算公式如下：

子带包络将用于PEAQ模型中的MOVs计算。

最后是响度计算部分，在响度计算中先是通过计算每帧音频信号的响度，最后求平均得到总响度，响度计算还包括噪声响度。

步骤1-5，模型输出参数(MOVs)的计算，包括根据上述四个步骤的计算结果(即参考信号和待测信号的激励样本、调制样本、频谱适应样本和子带包络、响度值、噪声值等参数)计算每帧音色感知特征参数值，最后通过平均计算得到待测信号的音色感知特征参数。

根据实验结果，本发明实施例使用的PEAQ模块提取的共有10个MOVs作为音色感知特征参数，即音色感知保真度的底层声学参数。这10个评价3D音频的音色感知特征参数有助于拟合主观测试结果，分别为：噪掩比TotalNMR、加窗后的调制差异WinModDiff、平均调制差异AvgModDiff、响度失真RmsNoiseLoud、参考信号带宽BandwidthRef、待测信号带宽Bandwidthtest、平均块失真ADB、差错信号谐波结构EHS、最大检测概率MFPD和相对干扰帧RelDisFrames。

对于每一个MOVs的详细计算过程为现有技术，可参考ITU的BS.1387-1标准。在本实施例中分别以平均调制差异AvgModDiff计算，待测信号带宽Bandwidthtest和参考信号带宽BandwidthRef为例说明MOVs的计算过程。

待测信号平均调制差异AvgModDiff，也称为暂态包络差异，也即是调制差异是通过计算每个巴克带滤波通道的调制差异也即是本地调制差异(local modulation differences)，并求平均可得，如下式所示：

$\mathrm{ModDiff}[k_1,n_1]=w\·\frac{\left|{\mathrm{Mod}}_{\mathrm{test}}\right[k_1,n_1]-{\mathrm{Mod}}_{\mathrm{Ref}}[k_1,n_1]}{\mathrm{offset}+{\mathrm{Mod}}_{\mathrm{Ref}}[k_1,n_1]}$

.........(10)

$\left|\begin{array}{cc}w=1.0& \left|{\mathrm{Mod}}_{\mathrm{test}}\right[k_1,n_1]>{\mathrm{Mod}}_{\mathrm{Ref}}[k_1,n_1]\\ w=0.1& \left|{\mathrm{Mod}}_{\mathrm{test}}\right[k_1,n_1]<{\mathrm{Mod}}_{\mathrm{Ref}}[k_1,n_1]\end{array}\right.$

其中，Mod_test[k₁，n₁]和Mod_Ref[k₁，n₁]是在步骤1-4具体计算中已经提到的待测信号和参考信号的子带包络，其中，offset表示偏移量，一般设定为0.01，w是加权因子，根据参考信号和待测信号子带包络相对大小而确定。那么对瞬态调制差异求平均即可得到所求平均调制差异：

$\mathrm{ModDiff}\left[n_1\right]=\frac{100}{Z_1}\underset{k_1=0}{\overset{Z_1-1}{\mathrm{\&Sigma;}}}\mathrm{ModDiff}[k_1,n_1]...\left(11\right)$

待测信号带宽Bandwidthtest、参考信号带宽BandwidthRef参数由参考信号和待测信号的FFT输出(以dB为单位)计算而来。具体方法为：

(1)在待测信号中，找出21.6kHz以上信号的最大幅值(单位dB)作为参考门限。

(2)在参考信号中，从21.6kHz开始往低频搜索，找到第一个超过参考门限10dB的样点，记下此样点的索引号K_R[n₁]，作为参考信号的带宽。

(3)在待测信号中，从上步记下的索引K_R[n₁]号开始往低频搜索，找到第一个超过参考门限5dB的样点，记下此时的索引号为K_R[n₁]，作为待测信号的带宽。

取K_R[n₁]的时域平均即为两个输入信号的最终带宽，上面计算的是瞬态的带宽，对所有帧求平均后，公式如下：

$W_T=\frac{1}{N_1}\underset{n=0}{\overset{N_1-1}{\mathrm{\&Sigma;}}}{K}_R\left[n_1\right]...\left(12\right)$

W_T表示信号整体平均带宽，N₁表示帧数，当然，这里要注意，如果通过上述方法计算的带宽没有超过8.1khz，则忽略该帧数据的带宽。

图3是本发明实施例中的水平声场感知模型框图，参考信号和待测信号的双耳合成信号s_test和s_ref输入环绕双耳模型，经过滤波和频带扩展输出人耳基底膜的激励模式；然后根据认知模型得到水平声场感知特征参数。

为便于实施参考起见，提供实施例计算水平声场感知特征参数的具体步骤如下所述：

步骤2-1，通过FFT变换将参考信号和待测信号的双耳合成信号变换为频域信号，将参考信号和待测信号的频域信号通过一个Gammatone滤波器组划分为24个频段的巴克带，模拟外耳和中耳滤波对划分的巴克带进行加权，得到参考信号的激励模式和待测信号的激励模式。

本发明实施例使用4阶的Patterson-Holdsworth滤波器组(也称为Gammatone滤波器组)，与PEAQ中的基于FFT的人耳感知模型相比，Gammatone滤波器组能很好地模拟听觉实验中的生理数据，能够体现基底膜尖锐的滤波特性。考虑到Gammatone滤波器组与人耳基底膜两者在冲激响应、幅频特性等方面的相似性，可以采用一组中心频率呈对数均匀分布的Gammatone滤波器组来模拟基底膜。

Gammatone滤波器用一个因果的冲激响应函数Gammatone函数来描述其滤波器特性，其时域表示如式13所示：

g(t)＝Bⁿt^n-1e^-2πBtcos(2πf_ot+φ)u(t).........(13)

其中：g(t)是Gammatone函数缩写，t是时间变量，φ表示Gammatone滤波器的初始相位。t＜0时u(t)＝0，t＞0时u(t)＝1，n为Gammatone滤波器的阶数，参数B＝b×ERB(f_o)，b₁＝1.019是为了让Gammatone函数更好地与生理数据相符而引入的参数，ERB(f_o)为Gammatone滤波器的等价矩形带宽，它同Gammatone滤波器中心频率f_o的关系是：

ERB(f_o)＝24.7+0.108f_o.........(14)

不同的特征频率对应不同的Gammatone滤波器中心频率，这个关系可以把不同Gammatone滤波器和不同的基底膜位置对应起来，从而实现对基底膜每个位置质点运动状态的仿真。使用Gammatone滤波器组对双耳合成信号进行滤波，则输出结果为待测信号的激励模式和参考信号的激励模式，可以分别用

表示，式中k₂和n₂分别表示Gammatone滤波器组频带和时间的索引。

步骤2-2，计算参考信号的激励模式和待测信号的激励模式的归一化互相关函数，得到参考信号和待测信号之间的双耳时间差异、参考信号和待测信号之间的双耳互相关系数，通过计算参考信号的激励模式和待测信号的激励模式的声强比值的对数值，得到参考信号和待测信号之间的双耳声级差异。

所述待测信号和参考信号的激励模式输入认知模型，得到表征水平声场感知特征参数：双耳声级差异(ILD)，双耳时间差异(ITD)，双耳互相关系数(IACC)。参考信号和待测信号的ITD主要用于低频(低于1500HZ)声源方向定位，ITD可以通过计算基于时间窗的归一化互相关函数(Normalized Cross-correlation Function：NCF)来计算，参考信号和待测信号的ITD，IACC的计算方法是一样的，以待测信号为例，互相关函数的计算是用长度约为20ms的7/8重叠矩形窗，待测信号的互相关函数如下式：

${\mathrm{NCF}}_{k_2,n_2}\left[d\right]=\frac{\underset{l}{\mathrm{\&Sigma;}}{X}_{L,k_2,n_2}\left[l\right]X_{R,k_2,n_2}[l+d]}{\sqrt{\underset{l}{\mathrm{\&Sigma;}}{X}_{L,k_2,n_2}^2\left[l\right]X_{R,k_2,n_2}^2\left[l\right]}}...\left(15\right)$

上式中d表示样本间的时间延时，l表示离散序列索引，

和

分别表示待测信号左右声道频谱序列，

表示待测信号的频谱序列经过d延时后的结果。

则IACC定义为所有计算样本间延时d的NCF最大值，而ITD即是该最大值所对应的d值，也即是：

$\mathrm{IACC}[k_2,n_2]=\max {\left|{\mathrm{NCF}}_{k_2,n_2}\right[d\left]\right|}_{d=-N_2}^{d=+N_2}...\left(16\right)$

$\mathrm{ITD}[k_2,n_2]\underset{d}{\arg }\max {\left|{\mathrm{NCF}}_{k_2,n_2}\right[d\left]\right|}_{d=-N_2}^{d=+N_2}...\left(17\right)$

其中，IACC[k₂，n₂]和ITD[k₂，n₂]分别表示双耳互相关系数序列和双耳时间差序列。N₂是样本延时d的范围，理论上几乎涵盖了所有ITD值，以样本数来表示。可感知的声源方向变化可以用基于单位圆的两个位置之间的欧式距离来近似计算。一般的，对于两个不同的方位角θ₁、θ₂之间的距离可以用以下公式计算：

$\sqrt{{({\cos \mathrm{\&theta;}}_1-{\cos \mathrm{\&theta;}}_2)}^2+{({\sin \mathrm{\&theta;}}_1-{\sin \mathrm{\&theta;}}_2)=}^2}\sqrt{2-2\cos ({\mathrm{\&theta;}}_1-{\mathrm{\&theta;}}_2)}...\left(18\right)$

从而，根据待测信号和参考信号之间的ITD的差别，两声源方向的感知距离计算如下：

$\mathrm{\&Delta;ITD}[k_2,n_2]=\sqrt{2-2\cos \mathrm{\&pi;}{f}_s/N_{\max }\left({\mathrm{ITD}}_{\mathrm{test}}\right[k_2,n_2]-{\mathrm{ITD}}_{\mathrm{ref}}[k_2,n_2\left]\right)}...\left(19\right)$

其中，ITD_test[k₂，n₂]和ITD_ref[k₂，n₂]表示待测信号和参考信号的双耳时间差序列，f_s是采样率，而N_max是样本间ITD的最大值，因此

可视作限制上述公式中余弦函数的输入在0-π之间的归一化因子。

在某些情况下，如果IACC值太小以至于不能产生稳定的可感知声源方向，ITD值便不易计算。因此，我们需要加入一个确定性因子p[k₂，n₂]来考虑可计算的ITD确定性，该确定性可通过经近似相切反曲函数处理的IACC非线性变换建模得到：

$p[k_2,n_2]={\{1+e^{-S\left(\right|\mathrm{IACC}[k_2,n_2]\left|\right)-T_k}\}}^{-1}...\left(20\right)$

S和T_k分别是斜度和阀值，T_k在不同的频带取不同的值。因此经过确定性因子的建模，ITD失真值ITDDist[k₂，n₂]计算如下：

$\mathrm{ITDDist}[k_2,n_2]=\frac{1}{2}\left(p_{\mathrm{test}}\right[k_2,n_2]+p_{\mathrm{ref}}[k_2,n_2\left]\right)\&times;\mathrm{\&Delta;ITD}[k_2,n_2]...\left(21\right)$

其中，p_ref[k₂，n₂]和p_test[k₂，n₂]分别表示参考信号和待测信号的确定性因子。

IACC表示清晰声源宽度，根据上述计算的IACC[k₂，n₂]，不难计算IACC失真值IACCDist[k₂，n₂]，与计算ITD不同这里计算互相关系数使用的窗函数长度为50ms：

IACCDist[k₂，n₂]＝w₂×|IACC_test[k₂，n₂]-IACC_ref[k₂，n₂]|.........(22)

其中，IACC_test[k₂，n₂]和IACC_ref[k₂，n₂]分别表示待测信号和参考信号的双耳互相关系数、w₂表示加权系数。

参考信号和待测信号的ILD主要用于高频(高于2500HZ)声源方向定位，用于描述左右耳接收到声压级的差异，ILD是高频声音方向感知测度的重要线索，因此，待测信号和参考信号之间的ILD差异可以反映高频音频信号可感知的声音方向降质程度。ILD的计算公式如下：

$\mathrm{ILD}[k_2,n_2]={10\log }_{10}\left(\frac{\underset{l}{\mathrm{\&Sigma;}}{X}_{L,k_2,n_2}^2\left[l\right]}{\underset{l}{\mathrm{\&Sigma;}}{X}_{R,k_2,n_2}^2\left[l\right]}\right)...\left(23\right)$

ILD值ILD[k₂，n₂]为10倍于第n₂时间帧第k₂频带的左右耳输入信号X_R和X_L强度比值的对数值，且强度值是由长度为10ms的3/4重叠矩形窗计算。ILD的计算区别于ITD，首先，根据ILD不同，两处声源的可感知差异可由ILD的对数域值(DB功率)来表征；其次，针对基于ILD的声源方向定位判断，会考虑声音强度的权重。因此ILD失真测度ILDDist[k₂，n₂]如下所示：

$\mathrm{ILDDist}[k_2,n_2]=w_3\left[k_2\right]\&times;{\log }_{10}\left(\underset{l}{\mathrm{\&Sigma;}}{X}_{T,k_2,n_2}^2\right[l]\&times;|{\mathrm{ILD}}_{\mathrm{test}}[k_2,n_2]-{\mathrm{ILD}}_{\mathrm{ref}}[k_2,n_2]\left|\right)...\left(24\right)$

其中，w₃[k₂]是非线性权重因子，表示映射每个频带的ILD失真相对重要度，ILD_test[k₂，n₂]和ILD_ref[k₂，n₂]分别表示待测信号和参考信号的双耳声级差。综合上述的ILD，ITD和IACC将作为水平声场感知特征参数输入到3D心理认知模型中处理。

图4是本发明实施例的高度感知模型结构框图，待测信号的双耳合成信号s_test和s_ref经过声道分离，分别以单声道(包括左右声道)作为高度感知模型的输入信号，然后经过FFT变换，高斯滤波，波峰和波谷检测，并采用估计器与HRTF数据库中的高度角信息进行对比，确定一次实验中待测信号的高度感知特征参数：声源的高度角θ。

为便于实施参考起见，提供实施例计算高度感知特征参数的具体步骤如下所述：

步骤3-1，分离待测信号的双耳合成信号的左右声道，并分别进行快速傅里叶变换，得到待测信号的频谱信号；

步骤3-2，将待测信号的频谱信号通过一个高斯滤波器滤波，得到滤波后的频谱信号；

步骤3-3，在待测信号的滤波后的频谱信号中检测频谱的波峰和波谷，将出现频谱的波峰波谷的位置和头相关传输函数数据库中的不同高度角的波峰波谷位置进行对比，估计得到待测信号的高度角。

在高度角θ的计算中，需要提供一个参考的数据库，数据库中包含头相关传输函数(HeadRelated Transfer Function，HRTF).估计器中用到的HRTF数据库是KEMAR数据库，该数据库来自麻省理工学院媒体实验室在假人工头(KEMAR)上测得的HRTF数据。在KEMAR模拟系统中搭建一个人工头放置在距离扬声器中垂线1.5米的位置接受双耳输入信号，人工头悬空并对称，左右耳输入信号假设是一致的，由估计器来估计声音事件的高度角。HRTFs是在无回声和无混响的听音室里，在中垂面上在-40°和90°之间以没10°为步长进行测量的，在KEMAR数据库有44100Hz的采样率和65个采样点的谱分辨率。

用于定位的谱线索主要集中于4k～16kHz，且决定前向声场方位的谱线索包含一个在4k～8kHz之间低于截止频率的八度音波谷(定义为N1)和一个能量增加且频率高于13kHz的八音度波谷；而上方声场方位的谱线索是在7k～9kHz之间的1/4八度音波峰(定义为P1)。随着声源的高度角在-45°和45°之间变化的时候，频谱的波谷在6k～10kHz之间也随着变化。多次实验的结果表明头传输函数的参数主要由第一个和第二个频谱波谷和第一个频谱波峰组成。

在频谱的波峰和波谷检测中，在5kHz以上的波峰和波谷对于感知高度角的贡献是比较大的，因此提取在频率范围4kHz以上的谱峰，根据具体位置对波峰和波谷进行标注(如P1，N1，P2，N2，....)。在提取波峰和波谷之前，应该用高斯滤波器对谱波动进行平滑处理：

$H_w\left(k_3\right)=\underset{n_3=-N_3}{\overset{N_3}{\mathrm{\&Sigma;}}}H(k_3+n_3)W\left(n_3\right)...\left(25\right)$

$W\left(n_3\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}\mathrm{\&sigma;}}{e}^{\frac{-{n_3}^2}{{2\mathrm{\&sigma;}}^2}}...\left(26\right)$

H(k₃+n₃)是不同高度角的频谱，H_w(k₃)是经过高斯滤波器后的频谱，W(n₃)是高斯滤波器，k₃和n₃表征离散频点，为了使得信号更平滑，式中的n₃和σ的设置依赖于采样率，在KEMAR曲线中N₃和σ分别赋值为4和1.3，所述的频谱的波峰和波谷被定义为水平谱的相对极大值和极小值，且其中水平谱被归一化至0dB以使得峰值和谷值和声音音量相互独立。

在高度角θ估计中，为了能从HRTF数据库中的抽取出和实际单耳信号谱匹配度最高的HRTF曲线，需要对比不同HRTF曲线的波峰和波谷从而在数据库里找到最匹配的高度角θ对应的HRTF曲线。估计器计算过程如下描述，每个单耳信号和每个HRTF谱都有一套波峰属性和波谷属性：N和P，而k代表频带索引，因此有

Signal＝{k_Ni，k_Pi}.........(28)

其中，

表示不同高度角的HRTF谱，k_Ni和k_Pi分别表示波峰和波谷索引，

表示不同高度角，Signal表示一组波峰和波谷。

上式中，由于体现了声源高度定位最重要的频率段是在4000Hz和15000Hz之间，因此只有在4000Hz和15000Hz之间的波峰和波谷用来分析估计。每个HRTF曲线都赋予一个

值，用来表示各自声源信号频谱的波峰和波谷的频谱距离，也即是

上式中，

为声源信号频谱的波峰和波谷的频谱距离，w_i表示加权因子。设k_i和k_j表示波峰或波谷的索引，

若k_i＝k_j则w_i＝1；若|k_i-k_j|＝1则w_i＝0.5；若|k_i-k_j|＝2则w_i＝0.25；其他情况下，w_i＝0。w值在两个波峰频率相等时取最大值，且随着两个波峰频率不同而下降。信号谱中的每个波峰和波谷都用来从HRTF中搜索最相近的波峰和波谷，在算法中计算HRTF数据库中每个高度角对应波谷和声源波谷的频谱距离

最后

的最大值所对应的HRTF高度角即是所估计的高度角θ，也即是

如果有多个角度有同样的

值，则在该时间窗下的高度角θ为这些角度的平均值。

步骤3，随机选取若干组步骤2所得待测信号的音色感知特征参数、水平声场感知特征参数、高度感知特征参数作为3D声场心理认知模型的输入，以预设的相应主观评价分数作为3D声场心理认知模型的输出，经过训练得到3D声场心理认知模型，将步骤2所得其他待测信号的音色感知特征参数、水平声场感知特征参数、高度感知特征参数输入训练后的3D声场心理认知模型，得到音频质量客观评价结果。随机选取的比例可以由本领域技术人员自行指定，一般是取一半待测信号的的音色感知特征参数、水平声场感知特征参数、高度感知特征参数作为3D声场心理认知模型的输入进行训练，将另一半待测信号的音色感知特征参数、水平声场感知特征参数、高度感知特征参数输入训练后的3D声场心理认知模型，得到音频质量客观评价结果。

图5是本发明实施例的3D声场心理认知模型结构框图，将之前所述的音色感知模型、水平声场感知模型和高度感知模型所得到的音色感知特征参数、水平声场感知特征参数、高度感知特征参数共14个参数作为BP神经网络的输入，经过神经网络后得到最终单一的音频质量客观评价分数。

实施例通过运用人工BP神经网络模型实现融合，基本实施过程如下：建立BP神经网络模型，其结构如图5所示，图中N表示节点。节点激活函数采用sigmoid函数，即

其中e为数学常数，x为该函数的输入。

该模型包含14个输入，并且神经网络模型的隐层中有5个节点。预先建立序列集，将原始的12个参考序列作为参考信号，利用不同的多声道压缩编码系统对12个序列进行编解码，得到12×12＝144个作为待测信号的测试序列。可以预先由专家对每个测试序列进行主观评价，得到主观评价分数，本发明中作为已知量输入即可。随机取其中72个测试序列作为训练用的待测信号，将根据这72个待测信号和参考信号的双耳合成信号所得音色感知特征参数、水平声场感知特征参数、高度感知特征参数作为输入，以预设的相应主观评价分数作为输出，对模型映射过程中输入层的两个限制因子a_min[i]和a_max[i]、输入层的加权系数w_x[i]、输出层的加权系数w_y[j]以及输出层的两个限制因子b_min和b_max进行训练，得到一系列系数集合，3D声场心理认知模型形成。剩余的72个测试序列作为评价用的待测信号，使用训练过的神经网络进行客观质量评价得到客观评分，即将根据剩余72个待测信号和参考信号的双耳合成信号所得音色感知特征参数、水平声场感知特征参数、高度感知特征参数作为输入，训练过的神经网络的输出即是音频质量客观评价分数。

将音色特征参数、水平声场感知特征参数和高度感知特征参数共14个参数作为神经网络模型的输入a[i]输入给神经网络模型，映射到失真索引(Distortion Index：DI)：

$\mathrm{DI}=w_y\left[J\right]+\underset{j=0}{\overset{J-1}{\mathrm{\&Sigma;}}}\left(w_y\right[j]\&times;\mathrm{sig}\left(w_x\right[I,j]+\underset{i=0}{\overset{I-1}{\mathrm{\&Sigma;}}}{w}_x[i,j]\&times;\frac{x\left[i\right]-a_{\min }\left[i\right]}{a_{\max }\left[i\right]-a_{\min }\left[i\right]}))...\left(31\right)$

式中，I＝14，J＝5，则w_y[J]表示j＝5时候的输出层加权系数，w_x[I，j]表示i＝14时的输入层加权系数，x[i]表示输入，w_x[i，j]表示输入层加权系数；i取0，1，K，I-1，j取0，1，K，J-1，利用失真索引DI，最终计算出音频质量客观评价分数ODG：

ODG＝b_min+(b_max-b_min)×sig(DI).........(32)

该分数即为本发明所提出的评价3D音频质量的最终得分。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (6)

1.一种3D音频质量客观评价方法，其特征在于，包含以下步骤：

步骤1，将多声道信号的参考序列下混得到参考信号的双耳合成信号，将多声道信号的测试序列下混得到待测信号的双耳合成信号；

步骤2，根据步骤1所得参考信号的双耳合成信号和待测信号的双耳合成信号，计算得到待测信号的音色感知特征参数、水平声场感知特征参数、高度感知特征参数； 

步骤3，随机选取若干组步骤2所得待测信号的音色感知特征参数、水平声场感知特征参数、高度感知特征参数作为3D声场心理认知模型的输入，以预设的相应主观评价分数作为3D声场心理认知模型的输出，经过训练得到3D声场心理认知模型，将步骤2所得其他待测信号的音色感知特征参数、水平声场感知特征参数、高度感知特征参数输入训练后的3D声场心理认知模型，得到音频质量客观评价结果。

2.根据权利要求1所述的一种3D音频质量客观评价方法，其特征在于：在步骤1中，采用头相关传输函数实现下混，得到双耳合成信号。

3.根据权利要求1或2所述的一种3D音频质量客观评价方法，其特征在于：采用ITU-R BS.1387-1标准的感知音频质量客观评价模块实现计算待测信号的音色感知特征参数，包括如下步骤，

步骤1-1，通过FFT变换将双耳合成信号变换为频域信号，对频域信号计算平方幅值得到能量信号，模拟人外耳和中耳听觉效应对能量信号进行滤波处理，得到滤波后的能量信号；

步骤1-2，将滤波后的能量信号进行临界带分组，得到临界带分组信号；所述临界带分组处理模拟了人耳基底膜对不同频带信号的响应差异；

步骤1-3，根据心理声学理论对临界带分组信号进行临界带扩展得到激励样本，并计算得到调制样本；

步骤1-4，对激励样本经过时域和频域平滑处理，然后通过电平和模式调整，得到频谱适应样本，并计算掩蔽阈值；

步骤1-5，根据上述四个步骤的计算结果计算每帧音色感知特征参数值，最后通过平均计算得到待测信号的音色感知特征参数。

4.根据权利要求1或2所述的一种3D音频质量客观评价方法，其特征在于：所述待测信号的水平声场感知特征参数包括双耳声级差异、双耳时间差异和双耳互相关系数，计算水平声场感知特征参数包括如下步骤，

步骤2-1，通过FFT变换将双耳合成信号变换为频域信号，将参考信号和待测信号的频域信号通过一个Gammatone滤波器组划分为24个频段的巴克带，模拟外耳和中耳滤波对划分的巴克带进行加权，得到参考信号的激励模式和待测信号的激励模式；

步骤2-2，计算参考信号的激励模式和待测信号的激励模式的归一化互相关函数，得到参考信号和待测信号之间的双耳时间差异、参考信号和待测信号之间的双耳互相关系数，通过计算参考信号的激励模式和待测信号的激励模式的声强比值的对数值，得到参考信号和待测信号之间的双耳声级差异。

5.根据权利要求1或2所述的一种3D音频质量客观评价方法，其特征在于：所述待测信号的高度感知特征参数包括高度角，计算高度感知特征参数包括如下步骤，

    步骤3-1，分离待测信号的双耳合成信号的左右声道，并分别进行快速傅里叶变换，得到待测信号的频谱信号；

步骤3-2，将待测信号的频谱信号通过一个高斯滤波器滤波，得到滤波后的频谱信号；

步骤3-3，在待测信号的滤波后的频谱信号中检测频谱的波峰和波谷，将出现频谱的波峰波谷的位置和头相关传输函数数据库中的不同高度角的波峰波谷位置进行对比，估计得到待测信号的高度角。

6.根据权利要求1或2所述的一种3D音频质量客观评价方法，其特征在于：所述3D声场心理认知模型使用BP神经网络，在3D声场心理认知模型中BP神经网络使用最速下降法。

Images