CN104408040A - 头相关函数三维数据压缩方法与系统 - Google Patents

Abstract

本发明实施例提供一种头相关脉冲响应三维数据压缩处理方法。本发明方法，包括：将HRIR三维数据重建为最小相位化的HRIR；将所述最小相位化的HRIR解卷为低阶FIR形式的仰角平面滤波器、高阶FIR形式的水平角平面滤波器；将所述高阶FIR形式的水平角平面滤波器建模为低阶IIR形式的水平角平面滤波器；所述处理器将所述低阶FIR形式的仰角平面滤波器以及所述IIR形式的水平角平面滤波器存储至本地缓存。本发明实施例克服了现有技术中HRIR三维数据压缩率低的问题。

Description

头相关函数三维数据压缩方法与系统

技术领域

本发明实施例涉及电通信领域，尤其涉及一种头相关脉冲响应三维数据压缩方法与系统。

背景技术

人类的听觉是三维感知，即不仅能感觉出声音的音调、音强、音色，还能分辨声源方向与距离。为了使不在现场的听众也能获得必要的“临场感”，从而发明了三维音频技术。三维音频技术可再现原始声场全部或部分信息，使听者产生“身临其境”的聆听体验。三维音频可通过双耳声学技术来实现，即采用头相关传递函数(Head Related Transfer function，以下简称：HRTF)对单通道声源信号进行滤波处理，将得到具有方向感的声音使用耳机或一对扬声器进行重放。双耳声学技术可应用于消费电子、移动终端、多媒体系统、虚拟现实、盲人听觉导航、飞行训练等领域。

HRTF的时域表示称为头相关脉冲响应(Head Related Impulse Responses，以下简称：HRIR)。HRIR描述了声音从声源传播到耳内的过程中，人的躯干、头部、耳部等对声音信号具有衍射、反射、遮挡等效应，这可以等效为滤波作用。实际中，通常采用经过精确实验测量得到的HRIR数据库，这种数据库是一个三维矩阵，该数据库存储了大量声源空间位置的HRIR。其中每个HRIR的冲激响应持续时间较长，等效FIR滤波器阶数很高，并且每个人的HRIR都不尽相同。因此，要存储多位测量者的HRIR数据，需要消耗巨大的存储空间。当基于嵌入式设备实现三维虚拟声时，其存储空间更为紧张，成为三维音频技术应用的瓶颈。为此，需要采用HRIR三维数据压缩技术，以难以察觉的听觉定位性能损失为代价，来大幅度减少HRIR数据库的存储空间。但现有技术中，采用HRIR三维数据压缩技术存在有压缩率低的问题。

发明内容

本发明实施例提供一种头相关脉冲响应三维数据压缩方法与系统，以克服现有技术中HRIR三维数据压缩率低的问题。

本发明实施例提供了一种头相关脉冲响应三维数据压缩方法，包括：

处理器将头相关脉冲响应HRIR三维数据重建为最小相位化的HRIR；

所述处理器将所述最小相位化的HRIR解卷为低阶有限冲击响应滤波器FIR形式的仰角平面滤波器、高阶FIR形式的水平角平面滤波器；

所述处理器将所述高阶FIR形式的水平角平面滤波器建模为低阶IIR形式的水平角平面滤波器；

所述处理器将所述低阶FIR形式的仰角平面滤波器以及无限冲击响应滤波器IIR形式的水平角平面滤波器存储至本地缓存。

进一步地，所述处理器将头相关脉冲响应HRIR三维数据重建为最小相位化的HRIR，包括：

选取基准信号，并将所述基准信号与所有HRIR三维数据相关，获取相关信号最大点的位置，所述位置即为每个HRIR对应于基准信号的时延；

移除所述每个HRIR的时延点，与所述基准信号进行对齐。

进一步地，所述将所述最小相位化的HRIR解卷为低阶FIR形式的仰角平面滤波器、高阶FIR形式的水平角平面滤波器，包括：

步骤一、采用公式

$c=\left[\begin{array}{c}{D}_1\\ D_2\\ .\\ .\\ .\\ D_I\end{array}\right]\×\left[\begin{array}{c}{h}_1\\ h_2\\ .\\ .\\ .\\ h_I\end{array}\right]---\left(1\right)$

提取所述仰角平面的HRIR公用系数CF，其中，所述c为CF，所述h₁为HRIR系数，所述I为HRIR系数的个数，所述D₁为方位相关系数DF组合成矩阵的形式，采用公式

其中，所述N为N为向量d₁的长度，所述d₁(0)为HRIR的DF系数d₁的第一个分量；

步骤二、通过J个所述CF得到矩阵A∈R^M×J，其中，所述R为矩阵空间，所述每个CF系数的维数，所述J为CF系数个数，M为CF系数的长度；

将所述矩阵A作为所述水平角平面的HRIR专用系数DF，采用公式

$c=\left[\begin{array}{c}{D}_1\\ D_2\\ .\\ .\\ .\\ D_I\end{array}\right]\×\left[\begin{array}{c}{h}_1\\ h_2\\ .\\ .\\ .\\ h_I\end{array}\right]---\left(3\right)$

提取每个水平角平面的CF；

步骤三、通过I个所述CF得到矩阵E∈R^N×I，其中，所述N为DF系数的维数，所述R为矩阵空间，所述I为DF个数；

步骤四、将所述矩阵E作为所述仰角平面的DF，采用公式(3)提取每个仰角平面的CF；

步骤五、通过J个所述CF得到矩阵A∈R^M×J；

步骤六、判断所述矩阵A与矩阵E是否收敛，若否，则重复步骤二至步骤五，若是，则通过矩阵A与矩阵E可得到其中，所述A为携带水平角位置信息的水平角CF矩阵，所述E为携带仰角位置信息的仰角CF矩阵，其中，所述为卷积运算符，所述h_ij:为HRIR三阶张量矩阵元素，所述i为矩阵行数，所述j为矩阵列数。

进一步地，所述将所述高阶FIR形式的水平角平面滤波器建模为低阶IIR形式的水平角平面滤波器，包括：

计算无限冲击响应滤波器IIR建模重构误差

$E_i\left(z\right)=H_i\left(z\right)-{\hat{H}}_i\left(z\right)=H_i\left(z\right)-\frac{A_i\left(z\right)}{B\left(z\right)}=\frac{H_i\left(z\right)B\left(z\right)-A_i\left(z\right)}{B\left(z\right)},i=\mathrm{1,2},...,I---\left(4\right)$

其中，所述为IIR滤波器的z域冲击响应函数，所述H_i(z)为FIR滤波器的z域冲击响应函数，所述A_i(z)为IIR滤波器的零点表达式，所述B(z)为IIR滤波器的极点表达式，所述i为第i个FIR滤波器，所述I为FIR滤波器个数；

假设第j次迭代的B(z)已知，当第j+1次迭代时，通过公式

$E_i^{(j+1)}\left(z\right)=\frac{H_i\left(z\right)B^{(j+1)}\left(z\right)-A_i^{(j+1)}\left(z\right)}{B^{\left(j\right)}\left(z\right)},i=\mathrm{1,2},...,I---\left(5\right)$

计算B^(j+1)(z)与A_i ^(j+1)(z)；

通过Z反变换将公式(5)转换为矩阵

$\left(\begin{array}{ccccc}{G}^{\left(j\right)}\left(z\right)& 0& ...& 0& {F_1}^{\left(j\right)}\left(z\right)\\ 0& G^{\left(j\right)}\left(z\right)& ...& 0& {F_2}^{\left(j\right)}\left(z\right)\\ ...& ...& ...& ...& ...\\ 0& 0& ...& G^{\left(j\right)}\left(z\right)& {F_I}^{\left(j\right)}\left(z\right)\end{array}\right)\left(\begin{array}{c}-a_1^{(j+1)}\\ -a_2^{(j+1)}\\ .\\ .\\ .\\ -a_I^{(j+1)}\\ b^{(j+1)}\end{array}\right)=\left(\begin{array}{c}-{f_1}^{\left(j\right)}\\ -{f_2}^{\left(j\right)}\\ .\\ .\\ .\\ -{f_I}^{\left(j\right)}\end{array}\right)---\left(6\right)$

通过所述矩阵(6)求得IIR系数b和a，其中，所述G^(j)(z)为B^(j)(z)的倒数，所述F₁ ^(j)(z)为H_i(z)与B^(j)(z)的比值，所述f₁ ^(j)为F₁ ^(j)(z)的反z变换。

进一步地，所述迭代B(z)的初始值b(0)通过公式

c＝(X^TX)^-1X^Th       (7)

求得，其中， $c=\left(\begin{array}{c}{a}_1^{\left(0\right)}\\ a_2^{\left(0\right)}\\ .\\ .\\ .\\ a_I^{\left(0\right)}\\ b^{\left(0\right)}\end{array}\right),h=\left(\begin{array}{c}{h}_1\\ h_2\\ .\\ .\\ .\\ h_I\end{array}\right),X=\left(\begin{array}{cccccc}\mathrm{\Δ}& 0& ...& 0& 0& H_1\\ 0& \mathrm{\Δ}& ...& 0& 0& H_2\\ ...& ...& ...& ...& ...& ...\\ 0& 0& ...& 0& ...& H_I\end{array}\right)$

$\mathrm{\Δ}=\left(\begin{array}{cccc}\mathrm{\δ}\left(0\right)& 0& ...& 0\\ \mathrm{\δ}\left(1\right)& \mathrm{\δ}\left(0\right)& ...& 0\\ ...& ...& ...& ...\\ \mathrm{\δ}(L-1)& \mathrm{\δ}(L-2)& ...& \mathrm{\δ}(L-P-1)\end{array}\right)H_i=\left(\begin{array}{cccc}0& 0& ...& 0\\ h_i\left(0\right)& 0& ...& 0\\ h_i\left(1\right)& h_i\left(0\right)& ...& 0\\ ...& ...& ...& ...\\ h_i(L-2)& h_i(L-3)& ...& h_i(L-Q-1)\end{array}\right),$

P＝Q＝10，所述c为IIR滤波器零点与极点多项式的系数矩阵，所述X为变换矩阵，所述Δ为单位脉冲函数组成矩阵，所述H为HRIR系数组成矩阵，所述T为转置符号，所述δ(0)为单位脉冲函数，所述L为FIR滤波器长度，所述P为IIR滤波器零点个数，所述h_i(0)为第i个HRIR系数的第一个分量，所述I为HRIR系数个数。

本发明实施例一种头相关脉冲响应三维数据压缩系统，包括：

最小相位重建模块，用于将HRIR三维数据重建为最小相位化的HRIR；

解卷模块，用于将所述最小相位化的HRIR解卷为低阶FIR形式的仰角平面滤波器、高阶FIR形式的水平角平面滤波器；

建模模块，用于将所述高阶FIR形式的水平角平面滤波器建模为低阶IIR形式的水平角平面滤波器；

存储模块，用于将所述低阶FIR形式的仰角平面滤波器以及所述IIR形式的水平角平面滤波器存储至本地缓存。

本发明实施例将HRIR三维数据重建为最小相位化的HRIR；将所述最小相位化的HRIR解卷为低阶FIR形式的仰角平面滤波器、高阶FIR形式的水平角平面滤波器；将所述高阶FIR形式的水平角平面滤波器建模为低阶IIR形式的水平角平面滤波器；所述处理器将所述低阶FIR形式的仰角平面滤波器以及所述IIR形式的水平角平面滤波器存储至本地缓存。本发明实施例克服了现有技术中HRIR三维数据压缩率低的问题。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明头相关脉冲响应三维数据压缩处理方法的流程图；

图2a为本发明CIPIC数据库中双耳原始HRIR示意图；

图2b为本发明经过最小相位重构后的HRIR示意图；

图3为本发明最小相位HRIR三维数据库分组示意图；

图4为本发明每一组HRIR三维矩阵的2D-CFD解卷结果示意图；

图5a为本发明仰角平面的CF系数提取示意图；

图5b为本发明水平角平面的CF系数提取示意图；

图6为本发明HRIR三维数据压缩处理系统结构示意图；

图7为本发明头相关坐标系示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

图1为本发明头相关脉冲响应三维数据压缩处理方法的流程图，如图1所示，本实施例的方法可以包括：

步骤101、处理器将头相关脉冲响应HRIR三维数据重建为最小相位化的HRIR；

步骤102、所述处理器将所述最小相位化的HRIR解卷为低阶有限冲击响应滤波器FIR形式的仰角平面滤波器、高阶FIR形式的水平角平面滤波器；

步骤103、所述处理器将所述高阶FIR形式的水平角平面滤波器建模为低阶IIR形式的水平角平面滤波器；

步骤104、所述处理器将所述低阶FIR形式的仰角平面滤波器以及无限冲击响应滤波器IIR形式的水平角平面滤波器存储至本地缓存。

进一步地，所述处理器将头相关脉冲响应HRIR三维数据重建为最小相位化的HRIR，包括：

选取基准信号，并将所述基准信号与所有HRIR三维数据相关，获取相关信号最大点的位置，所述位置即为每个HRIR对应于基准信号的时延；

移除所述每个HRIR的时延点，与所述基准信号进行对齐。

进一步地，所述将所述最小相位化的HRIR解卷为低阶FIR形式的仰角平面滤波器、高阶FIR形式的水平角平面滤波器，包括：

步骤一、采用公式

$c=\left[\begin{array}{c}{D}_1\\ D_2\\ .\\ .\\ .\\ D_I\end{array}\right]\×\left[\begin{array}{c}{h}_1\\ h_2\\ .\\ .\\ .\\ h_I\end{array}\right]---\left(1\right)$

提取所述仰角平面的HRIR公用系数CF，其中，所述c为CF，所述h₁为HRIR系数，所述I为HRIR系数的个数，所述D₁为方位相关系数DF组合成矩阵的形式，采用公式

其中，所述N为N为向量d₁的长度，所述d₁(0)为HRIR的DF系数d₁的第一个分量；

步骤二、通过J个所述CF得到矩阵A∈R^M×J，其中，所述R为矩阵空间，所述每个CF系数的维数，所述J为CF系数个数，M为CF系数的长度；

将所述矩阵A作为所述水平角平面的HRIR专用系数DF，采用公式

$c=\left[\begin{array}{c}{D}_1\\ D_2\\ .\\ .\\ .\\ D_I\end{array}\right]\×\left[\begin{array}{c}{h}_1\\ h_2\\ .\\ .\\ .\\ h_I\end{array}\right]---\left(3\right)$

提取每个水平角平面的CF系数；

步骤三、通过I个所述CF得到矩阵E∈R^N×I，其中，所述N为DF系数的维数，所述R为矩阵空间，所述I为DF个数；

步骤四、将所述矩阵E作为所述仰角平面的DF，采用公式(3)提取每个仰角平面的CF系数；

步骤五、通过J个所述CF得到矩阵A∈R^M×J；

步骤六、判断所述矩阵A与矩阵E是否收敛，若否，则重复步骤二至步骤五，若是，则通过矩阵A与矩阵E可得到其中，所述A为携带水平角位置信息的水平角CF矩阵，所述E为携带仰角位置信息的仰角CF矩阵，所述为卷积运算符，所述h_ij:为HRIR三阶张量矩阵元素，所述i为矩阵行数，所述j为矩阵列数。

进一步地，所述将所述高阶FIR形式的水平角平面滤波器建模为低阶IIR形式的水平角平面滤波器，包括：

计算无限冲击响应滤波器IIR建模重构误差

$E_i\left(z\right)=H_i\left(z\right)-{\hat{H}}_i\left(z\right)=H_i\left(z\right)-\frac{A_i\left(z\right)}{B\left(z\right)}=\frac{H_i\left(z\right)B\left(z\right)-A_i\left(z\right)}{B\left(z\right)},i=\mathrm{1,2},...,I---\left(4\right)$

其中，所述为IIR滤波器的z域冲击响应函数，所述H_i(z)为FIR滤波器的z域冲击响应函数，所述A_i(z)为IIR滤波器的零点表达式，所述B(z)为IIR滤波器的极点表达式，所述i为第i个FIR滤波器，所述I为FIR滤波器个数；

假设第j次迭代的B(z)已知，当第j+1次迭代时，通过公式

$E_i^{(j+1)}\left(z\right)=\frac{H_i\left(z\right)B^{(j+1)}\left(z\right)-A_i^{(j+1)}\left(z\right)}{B^{\left(j\right)}\left(z\right)},i=\mathrm{1,2},...,I---\left(5\right)$

计算B^(j+1)(z)与A_i ^(j+1)(z)；

通过Z反变换将公式(5)转换为矩阵

$\left(\begin{array}{ccccc}{G}^{\left(j\right)}\left(z\right)& 0& ...& 0& {F_1}^{\left(j\right)}\left(z\right)\\ 0& G^{\left(j\right)}\left(z\right)& ...& 0& {F_2}^{\left(j\right)}\left(z\right)\\ ...& ...& ...& ...& ...\\ 0& 0& ...& G^{\left(j\right)}\left(z\right)& {F_I}^{\left(j\right)}\left(z\right)\end{array}\right)\left(\begin{array}{c}-a_1^{(j+1)}\\ -a_2^{(j+1)}\\ .\\ .\\ .\\ -a_I^{(j+1)}\\ b^{(j+1)}\end{array}\right)=\left(\begin{array}{c}-{f_1}^{\left(j\right)}\\ -{f_2}^{\left(j\right)}\\ .\\ .\\ .\\ -{f_I}^{\left(j\right)}\end{array}\right)---\left(6\right)$

通过所述矩阵(5)求得IIR系数b和a，其中，所述G^(j)(z)为B^(j)(z)的倒数，

所述F₁ ^(j)(z)为H_i(z)与B^(j)(z)的比值，所述f₁ ^(j)为F₁ ^(j)(z)的反z变换。

进一步地，所述迭代B(z)的初始值b(0)通过公式

c＝(X^TX)^-1X^Th          (7)

求得，其中， $c=\left(\begin{array}{c}{a}_1^{\left(0\right)}\\ a_2^{\left(0\right)}\\ .\\ .\\ .\\ a_I^{\left(0\right)}\\ b^{\left(0\right)}\end{array}\right),h=\left(\begin{array}{c}{h}_1\\ h_2\\ .\\ .\\ .\\ h_I\end{array}\right),X=\left(\begin{array}{cccccc}\mathrm{\Δ}& 0& ...& 0& 0& H_1\\ 0& \mathrm{\Δ}& ...& 0& 0& H_2\\ ...& ...& ...& ...& ...& ...\\ 0& 0& ...& 0& ...& H_I\end{array}\right)$

$\mathrm{\Δ}=\left(\begin{array}{cccc}\mathrm{\δ}\left(0\right)& 0& ...& 0\\ \mathrm{\δ}\left(1\right)& \mathrm{\δ}\left(0\right)& ...& 0\\ ...& ...& ...& ...\\ \mathrm{\δ}(L-1)& \mathrm{\δ}(L-2)& ...& \mathrm{\δ}(L-P-1)\end{array}\right)H_i=\left(\begin{array}{cccc}0& 0& ...& 0\\ h_i\left(0\right)& 0& ...& 0\\ h_i\left(1\right)& h_i\left(0\right)& ...& 0\\ ...& ...& ...& ...\\ h_i(L-2)& h_i(L-3)& ...& h_i(L-Q-1)\end{array}\right),$

P＝Q＝10，所述c为IIR滤波器零点与极点多项式的系数矩阵，所述X为变换矩阵，所述Δ为单位脉冲函数组成矩阵，所述H为HRIR系数组成矩阵，所述T为转置符号，所述δ(0)为单位脉冲函数，所述L为FIR滤波器长度，所述P为IIR滤波器零点个数，所述h_i(0)为第i个HRIR系数的第一个0分量，所述I为HRIR系数个数。

具体来说，最小相位重建的具体步骤为：

将每一个头相关脉冲响应(Head Related Impulse Responses，以下简称：HRIR)重建为最小相位化的HRIR，即去除每个HRIR的时延，也就是去除每个HRIR最开始幅值几乎为零的样本点，如图2a和图2b所示。实现该模块的方法根5据具体情况而定，当HRIR数据库中给出了声音到达双耳的时间，即各个HRIR的时延，那么就根据这组数据将各个HRIR的时延移除。当HRIR数据库没有给出各个HRIR的时延，那么就采用相关法来进行最小相位化重建。

首先，选取一个基准信号，与所有HRIR作相关。相关信号最大点的位置就是各个HRIR相对于基准信号的时延。将各个HRIR前面的时延点移除，进行对齐。基准信号的选取应选择与HRIR相关性强的，例如，将方位θ＝0°,φ＝0°处HRIR的前面几乎为零的时延点进行人工移除，然后将其作为基准信号。

公用因子分解CFD方法的目的为，对于一组给定的HRIR，表示为h_i，i＝1,2,…,I，对这组h_i进行解卷积，得到c与d_i，i＝1,2,…,I。其中，c称为CF系数，为这组HRIR公用系数；d_i称为DF系数，为每个HRIR的系数，与每个HRIR的方位相关。两者均为向量，且长度分别为M和N。h_i与c、d_i的关系为：

$h_i=c\⊗d_i,i=\mathrm{1,2},...,I---\left(8\right)$

其中，符号表示卷积。将式(8)写为矩阵相乘的形式，则有

h_i＝Cd_i，i＝1,2,…,I         (9)

其中，

其中，L为HRIR冲激响应长度，N为向量d_i的长度，M为向量c长度，三者之间的关系为L＝M+N-1。

根据式(9)，计算d(n)，则有

其中，表示表示矩阵的伪逆。

在得到d_i后，用下式更新c

$c=\left[\begin{array}{c}{D}_1\\ D_2\\ .\\ .\\ .\\ D_I\end{array}\right]\×\left[\begin{array}{c}{h}_1\\ h_2\\ .\\ .\\ .\\ h_I\end{array}\right]---\left(11\right)$

其中，D_i由d_i按照C的方式组合而成。

每当计算得到c后，将其带入式即可得到新的d_i，重复此过程，最终这组HRIR可解卷为CF与DF两部分系数，即c与d_i，i＝1,2,…,I。

对于最小相位化的HRIR三维矩阵对其进行分组，如图3所示，为将其平均分为2组的示意图。本发明将该三维矩阵按照此分组方法分为6组，每一组HRIR仍然是一个三维矩阵，将每一组HRIR水平方向上的HRIR个数记为J，仰角方向上HRIR个数记为I。分组后，对每一组HRIR进行2D-CFD分解，经过2D-CFD后，如图4所示，在每个仰角方向上得到一个水平角平面滤波器其长度为M，在每个仰角方向上得到一个仰角平面滤波器其长度为N。

处理器将最小相位化的HRIR解卷为低阶FIR形式的仰角平面滤波器、高阶FIR形式的水平角平面滤波器的(2D-CFD算法)描述如下：

(1)如图5a所示，根据式(1)提取每个仰角平面的CF系数，j＝1,2,…,J，DF系数初始值设为全1。将得到的J个组成矩阵A∈R^M×J；

(2)如图5b所示，将上一步得到的A作为每个水平角平面DF因子。采用式(2)，提取每个水平角平面CF系数i＝1,2，…,I。将得到的I个组成矩阵E∈R^N×I；

(3)如图5a所示，将上一步得到的E作为每个仰角平面的DF因子。采用式(4)，提取每个仰角平面的CF系数j＝1,2,…,J。将得到的J个CF组成矩阵A∈R^M×J；

(4)判断矩阵A与E是否收敛，若不收敛，则返回(2)-(3)，继续迭代；若收敛，算法结束。

其中，本实施例将M取为189，将N取为12。

经过上述迭代过程后，每一组HRIR三维数据集都被分解为两个矩阵，A∈R^M×J与E∈R^N×I。计算出A与E后，每一个HRIR可通过下式重建；

$h_{\mathrm{ij}:}=A(:,j)\⊗E(:,j)---\left(12\right)$

其中，表示卷积，A(:,i)为CF矩阵,E(:,J)为DF矩阵。

为进一步减少存储量，将FIR滤波器形式的HRIR建模为共极点IIR滤波器，即零极点模型形式CAPZ的滤波器。FIR滤波器的数学表达式为

$H_i\left(z\right)=\underset{l=0}{\overset{L-1}{\mathrm{\Σ}}}{c}_i\left(l\right)z^{-l}---\left(13\right)$

其中，L为滤波器阶数,l为第l阶自变量，z为z域自变量，c_i(l)为滤波器系数，

将其建模为IIR滤波器，数学表达式为

${\hat{H}}_i\left(z\right)=\frac{A_i\left(z\right)}{B\left(z\right)}=\frac{\underset{p=0}{\overset{P}{\mathrm{\Σ}}}{a}_i\left(p\right)z^{-p}}{1-\underset{q=1}{\overset{Q}{\mathrm{\Σ}}}b\left(q\right)z^{-q}}---\left(14\right)$

其中，P和Q分别为建模后的IIR滤波器的分子多项式阶数和分母多项式阶数，q为第q极点下标b(q)为第q个极点系数，a_i(p)为第p个零点系数。

将得到的6组水平角平面滤波器E1、E2、…、E6组，采用基于迭代线性最小二乘的方法。

图6为本发明HRIR三维数据压缩处理系统结构示意图，如图6所示，本实施例的系统，包括3个模块：最小相位重建模块、解卷模块、建模模块。

最小相位重建模块101，用于将HRIR三维数据重建为最小相位化的HRIR；

解卷模块102，用于将所述最小相位化的HRIR解卷为低阶FIR形式的仰角平面滤波器、高阶FIR形式的水平角平面滤波器；

建模模块103，用于将所述高阶FIR形式的水平角平面滤波器建模为低阶IIR形式的水平角平面滤波器；

存储模块104，用于将所述低阶FIR形式的仰角平面滤波器以及所述IIR形式的水平角平面滤波器存储至本地缓存。

具体来说，最小相位重建模块用于将HRIR重建为最小相位形式的信号，得到最小相位化的HRIR三维数据库；2D-CFD用于将三维HRIR数据集解卷为两组滤波器，一组为低阶FIR形式的仰角平面滤波器，另一组为高阶FIR形式的水平角平面滤波器；CAPZ模块用于将高阶FIR形式的水平角平面滤波器，建模为低阶IIR形式的水平角平面滤波器。已有的HRTF压缩方法压缩效果不够好，得到的数据存储量依然很大。本发明针对HRIR数据库耗费存储量大的问题，采用2D-CFD与基于迭代最小二乘的共极点IIR建模方法对HRIR数据库进行压缩，在相同的重构误差下，采用本发明对HRIR数据库压缩后所需要的存储量更小。

本发明实施例将HRIR三维数据重建为最小相位化的HRIR；将所述最小相位化的HRIR解卷为低阶FIR形式的仰角平面滤波器、高阶FIR形式的水平角平面滤波器；将所述高阶FIR形式的水平角平面滤波器建模为低阶IIR形式的水平角平面滤波器；所述处理器将所述低阶FIR形式的仰角平面滤波器以及所述IIR形式的水平角平面滤波器存储至本地缓存。本发明实施例克服了现有技术中HRIR三维数据压缩率低的问题。

为了验证本发明的有效性，进行了计算机仿真实验。实验时采用的数据库为加州大学戴维斯分校影像处理计算中心(Center for Image Processing andIntegrated Computing，以下简称：CIPIC)HRIR数据库中测量个体subject_003的双耳HRIR数据，并将双耳数据合并为一个三维矩阵，三维HRIR矩阵尺寸为D₁×D₂×L，数值为50×50×200。

(1)本方法的压缩率与重构误差

本方案的压缩率DR_2DCFD+CAPZ计算如下式所示，DR_2DCFD+CAPZ越小，所需要存储量越小，压缩效果越好。

${\mathrm{DR}}_{2\mathrm{DCFD}+\mathrm{CAPZ}}=\frac{[(P+1)\×D_1+Q]\×\mathrm{Nseg}+N\×D_2}{D_1\×D_2\×L}---\left(15\right)$

其中，P为IIR建模的零点个数，Q为IIR建模的极点个数，D₁为三维HRIR矩阵的第一维大小，N为DF系数的向量长度，D₂为三维HRIR矩阵的第二维大小，Nseg为分组数，本方案中将其设为6。

现有技术一的压缩率为

${\mathrm{DR}}_{1\mathrm{DCFD}+\mathrm{BMT}}=\frac{D_2\×M\×D_1+(P+1+Q)\×D_1}{D_1\×D_2\×L}---\left(16\right)$

现有技术二的压缩率与本方案相同。

文献己经证明，当重构误差小于10％时，HRIR可以保持其声学定位性能不变。据此，当重构误差最大但不超过10％、并且三种重构误差相同的情况下下，本发明与现有技术一、现有技术二的压缩率比较如表1所示。从表1可以看出，与现有技术一、现有技术二相比，本方法压缩率最低。表1为本发明与现有技术压缩率比较表。

表1

方案	压缩率
		1D-CFD+BMT(现有技术一)	3.21％
2D-CFD+基于最小二乘共极点IIR建模(现有技术二)	0.70％
		2D-CFD+基于迭代最小二乘共极点IIR建模(本发明)	0.60％

(2)主观听音测试

测试HRIR集：两组，分别为CIPIC的subject_003的HRIR集数据，及其压缩重构数据集。

测试主体：4名听觉正常的测试者：2名女性、2名男性；

测试角度(φ，θ)：对7个方位进行测试，分别为(0°，-35°)、(0°，0°)、(0°，35°)、(-22.5°，0°)、(33.8°，0°)、(61.9°，0°)、(90°，0°)。由于使用通用HRIR存在前后混淆这种可能性，因此选取的这7个方位均位于前半平面，头相关坐标系如图7所示。

测试音频：将采样率为44.1kHz的0.4秒白噪声与方位相同的两组HRIR分别卷积(原始HRIR、压缩后重建HRIR)，得到0.4秒的双耳声。将其重复10次，以0.1s的静音做间隔，得到该方位两组5秒测试音频。

将每个方位两组测试音频分别播放给4名测试者，进行听音比较，选择两组间角度差异。然后根据表2中差异角度查找对应得分。将7个方位的得分做平均，得到表3所示的听音测试结果。由表3可见，三种方案经过压缩、重构后的HRIR在听音方面与原始HRIR区别。表2为感知角度距离选择表。表3为听音测试结果表。

表2

差异角度	描述	相似分数
			0°	完全相同	5
0°～5°	刚好可察觉的差异	3
			5°～15°	显著的差异	2
15°～30°	非常不同	1
			>30°	完全不同	0

表3

(3)计算量

当需要某个方位的HRIR时，需要将最终压缩得到的水平方位IIR滤波器与仰角方位FIR滤波器进行卷积，重构每个方位IIR形式的HRIR，三个方案的计算量如表4所示。从表4可以看出，现有技术一的计算量最少，其次是本方案。表4为计算复杂度比较表。

表4

(1)解卷模块中，也可以采用其他的三维矩阵分解方式。

(2)建模模块中，也可以采用其他的共极点IIR建模方法，如基于均衡模型截断(Balanced Model Truncation，简称：BMT)的共极点IIR建模方法。

本发明采用2D-CFD方法，提高了HRIR三维数据集的压缩效果，降低了数据集的存储量。本发明采用迭代最小二乘方法进行共极点IIR建模，提高了HRIR三维数据集的压缩效果，降低了数据集的存储量，并减少了重构运算量。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (6)

1.一种头相关脉冲响应三维数据压缩方法，其特征在于，包括：

处理器将头相关脉冲响应HRIR三维数据重建为最小相位化的HRIR；

所述处理器将所述最小相位化的HRIR解卷为低阶有限冲击响应滤波器FIR形式的仰角平面滤波器、高阶FIR形式的水平角平面滤波器；

所述处理器将所述高阶FIR形式的水平角平面滤波器建模为低阶IIR形式的水平角平面滤波器；

所述处理器将所述低阶FIR形式的仰角平面滤波器以及无限冲击响应滤波器IIR形式的水平角平面滤波器存储至本地缓存。

2.根据权利要求1所述的方法，其特征在于，所述处理器将头相关脉冲响应HRIR三维数据重建为最小相位化的HRIR，包括：

选取基准信号，并将所述基准信号与所有HRIR三维数据相关，获取相关信号最大点的位置，所述位置即为每个HRIR对应于基准信号的时延；

移除所述每个HRIR的时延点，与所述基准信号进行对齐。

3.根据权利要求1所述的方法，其特征在于，所述将所述最小相位化的HRIR解卷为低阶FIR形式的仰角平面滤波器、高阶FIR形式的水平角平面滤波器，包括：

步骤一、采用公式

$c=\left[\begin{array}{c}{D}_1\\ D_2\\ .\\ .\\ .\\ D_I\end{array}\right]\×\left[\begin{array}{c}{h}_1\\ h_2\\ .\\ .\\ .\\ h_I\end{array}\right]---\left(1\right)$

提取所述仰角平面的HRIR公用系数CF，其中，所述c为CF，所述h₁为HRIR系数，所述I为HRIR系数的个数，所述D₁为方位相关系数DF组合成矩阵的形式，采用公式

其中，所述N为向量d₁的长度，所述d₁(0)为HRIR的DF系数d₁的第一个分量；

步骤二、通过J个所述CF得到矩阵A∈R^M×J，其中，所述R为矩阵空间，所述每个CF系数的维数，所述J为CF系数个数，M为CF系数的向量长度；

将所述矩阵A作为所述水平角平面的HRIR专用系数DF，采用公式

$c=\left[\begin{array}{c}{D}_1\\ D_2\\ .\\ .\\ .\\ D_I\end{array}\right]\×\left[\begin{array}{c}{h}_1\\ h_2\\ .\\ .\\ .\\ h_I\end{array}\right]---\left(3\right)$

提取每个水平角平面的CF；

步骤三、通过I个所述CF得到矩阵E∈R^N×I，其中，所述N为DF的维数，所述R为矩阵空间，所述I为DF个数；

步骤四、将所述矩阵E作为所述仰角平面的DF，采用公式(3)提取每个仰角平面的CF；

步骤五、通过J个所述CF得到矩阵A∈R^M×J；

步骤六、判断所述矩阵A与矩阵E是否收敛，若否，则重复步骤二至步骤五，若是，则通过矩阵A与矩阵E可得到其中，所述A为携带水平角位置信息的水平角CF矩阵，所述E为携带仰角位置信息的仰角CF矩阵，所述为卷积运算符，所述h_ij:为HRIR三阶张量矩阵元素，所述i为矩阵行数，所述j为矩阵列数。

4.根据权利要求1所述的方法，其特征在于，所述将所述高阶FIR形式的水平角平面滤波器建模为低阶IIR形式的水平角平面滤波器，包括：

计算无限冲击响应滤波器IIR建模重构误差

$E_i\left(z\right)=H_i\left(z\right)-{\hat{H}}_i\left(z\right)=H_i\left(z\right)-\frac{A_i\left(z\right)}{B\left(z\right)}=\frac{H_i\left(z\right)B\left(z\right)-A_i\left(z\right)}{B\left(z\right)},i=\mathrm{1,2},...,I---\left(4\right)$

其中，所述为IIR滤波器的z域冲击响应函数，所述H_i(z)为FIR滤波器的z域冲击响应函数，所述A_i(z)为IIR滤波器的零点表达式，所述B(z)为IIR滤波器的极点表达式，所述i为第i个FIR滤波器，所述I为FIR滤波器个数；

假设第j次迭代的B(z)已知，当第j+1次迭代时，通过公式

$E_i^{(j+1)}\left(z\right)=\frac{H_i\left(z\right)B^{(j+1)}\left(z\right)-A_i^{(j+1)}\left(z\right)}{B^{\left(j\right)}\left(z\right)},i=\mathrm{1,2},...,I---\left(5\right)$

计算B^(j+1)(z)与A_i ^(j+1)(z)；

通过Z反变换将公式(5)转换为矩阵

$\left(\begin{array}{ccccc}{G}^{\left(j\right)}\left(z\right)& 0& ...& 0& {F_1}^{\left(j\right)}\left(z\right)\\ 0& G^{\left(j\right)}\left(z\right)& ...& 0& {F_2}^{\left(j\right)}\left(z\right)\\ ...& ...& ...& ...& ...\\ 0& 0& ...& G^{\left(j\right)}\left(z\right)& {F_I}^{\left(j\right)}\left(z\right)\end{array}\right)\left(\begin{array}{c}-a_1^{(j+1)}\\ -a_2^{(j+1)}\\ .\\ .\\ .\\ -a_I^{(j+1)}\\ b^{(j+1)}\end{array}\right)=\left(\begin{array}{c}-{f_1}^{\left(j\right)}\\ -{f_2}^{\left(j\right)}\\ .\\ .\\ .\\ -{f_I}^{\left(j\right)}\end{array}\right)---\left(6\right)$

通过所述矩阵(6)求得IIR系数b和a，其中，所述G^(j)(z)为B^(j)(z)的倒数，所述F₁ ^(j)(z)为H_i(z)与B^(j)(z)的比值，所述f₁ ^(j)为F₁ ^(j)(z)的反z变换。

5.根据权利要求4所述的方法，其特征在于，所述迭代B(z)的初始值b(0)通过公式

c＝(X^TX)^-1X^Th         (7)

求得，其中， $c=\left(\begin{array}{c}{a}_1^{\left(0\right)}\\ a_2^{\left(0\right)}\\ .\\ .\\ .\\ a_I^{\left(0\right)}\\ b^{\left(0\right)}\end{array}\right),h=\left(\begin{array}{c}{h}_1\\ h_2\\ .\\ .\\ .\\ h_I\end{array}\right),X=\left(\begin{array}{cccccc}\mathrm{\Δ}& 0& ...& 0& 0& H_1\\ 0& \mathrm{\Δ}& ...& 0& 0& H_2\\ ...& ...& ...& ...& ...& ...\\ 0& 0& ...& 0& ...& H_I\end{array}\right)$

$\mathrm{\Δ}=\left(\begin{array}{cccc}\mathrm{\δ}\left(0\right)& 0& ...& 0\\ \mathrm{\δ}\left(1\right)& \mathrm{\δ}\left(0\right)& ...& 0\\ ...& ...& ...& ...\\ \mathrm{\δ}(L-1)& \mathrm{\δ}(L-2)& ...& \mathrm{\δ}(L-P-1)\end{array}\right)H_i=\left(\begin{array}{cccc}0& 0& ...& 0\\ h_i\left(0\right)& 0& ...& 0\\ h_i\left(1\right)& h_i\left(0\right)& ...& 0\\ ...& ...& ...& ...\\ h_i(L-2)& h_i(L-3)& ...& h_i(L-Q-1)\end{array}\right),$

P＝Q＝10，所述c为IIR滤波器零点与极点多项式的系数矩阵，所述X为变换矩阵，所述Δ为单位脉冲函数组成矩阵，所述H为HRIR系数组成矩阵，所述T为转置符号，所述δ(0)为单位脉冲函数，所述L为FIR滤波器长度，所述P为IIR滤波器零点个数，所述h_i(0)为第i个HRIR系数的第一个分量，所述I为HRIR系数个数。

6.一种头相关脉冲响应三维数据压缩系统，其特征在于，包括：

最小相位重建模块，用于将HRIR三维数据重建为最小相位化的HRIR；

解卷模块，用于将所述最小相位化的HRIR解卷为低阶FIR形式的仰角平面滤波器、高阶FIR形式的水平角平面滤波器；

建模模块，用于将所述高阶FIR形式的水平角平面滤波器建模为低阶IIR形式的水平角平面滤波器；

存储模块，用于将所述低阶FIR形式的仰角平面滤波器以及所述IIR形式的水平角平面滤波器存储至本地缓存。