CN102789500B - 一种音频比较方法 - Google Patents

Abstract

本发明公开了一种快速音频比较的方法，其所述音频比较部分包括以下步骤：读取音频p和音频q，划分出音频p的特征段集合Cp和音频q的特征段集合Cq；利用GPU(显示卡)上CUDA(通用设备计算架构)功能快速的计算出音频p的特征段集合Cp的能量特征值序列Wp和音频q的特征段集合Cq的能量特征值序列Wq；将各个特征段的特征值按照特征段的顺序构成能量矩阵；利用后缀数组的变形算法找出两个特征值序列的公共特征段集合Seg；快速扫描公共特征段集合Seg，找出连通区域，返回连通区域的集合Vres作为音频比较的结果；在波形图上标识比较结果。

Description

一种音频比较方法

技术领域

本发明涉及一种音频比较方法，属于计算机音频、多媒体信息技术以及模式匹配等领域，特别涉及一种在音频数据库中基于能量信息提取特征的两段音频数据的快速比较。

背景技术

随着计算机技术的发展，特别是海量存储设备和大容量内存在PC机上的实现，对音频媒体进行数字化处理便成为可能。数字化处理的核心是对音频信息的采样，通过对采集到的样本进行加工，达成各种效果，这是音频媒体数字化处理的基本含义。为了检查处理过后的音频文件和原始文件之间有哪些区别，如何快速、准确地从两段音频中找到相同和不同之处是影视后期编辑中最基本也是最重要的问题之一。解决这个问题的有效手段便是音频比较技术。常用的音频技术主要可分为基于采样数据的和基于音频内容的二种。其中，基于采样数据的音频比较技术通常通过直接在比较两段音频的每个采样点来实现。这种技术存在两个明显不足：首先，音频采样点数据的海量性决定了对于两段较长的音频文件逐点比较这样的问题是NP的，不可能很快的得出结果，没有现实性的意义。其次，基于逐点比较的方法会使音频中存在大量的破碎段落，使得比较的结果过于零碎，不能提供直观可见的结果。基于音频内容的音频比较技术主要是对音频本身的各种特征进行分析和比较的技术。目前常用的音频特征有：波形、音质和MFCC等。但基于波形的比较方法过于依赖波形图的拟合方程，使得比较结果的存在很大误差。而基于的音质方法存在音质标准和接受器材约束的问题，使得比较的精度和效率也都不够好。基于MFCC特征的方法由于需要将音频按照语句划分出来，无法大规模快速进行。还有些音频比较方法采用了半监督的反馈技术，通过与用户的交互来改善比较的结果，直到用户满意。这种方法通常会给用户带来较大的负担和时间开销。

发明内容

发明目的：本发明所要解决的技术问题是针对现有技术的不足，并且为影视后期制作公司提供一个新的解决方案，提供一种音频比较方法。

技术方案：本发明公开了一种音频方法，所述音频快速比较部分包括以下步骤：

步骤1，将待比较的音频p和音频q，分别按照其音频的长度划分出长度相等的两组音频数据特征段，其中音频p划分为特征段集合Cp，Cp={cp₁，cp₂…cp_i…cp_Lp}，音频q划分为特征段集合Cq，Cq={cq₁，cq₂…cq_j…cq_Lq}，其中cp_i表示特征段集合Cp中第i个音频数据特征段，1≤i≤Lp，Lp为特征段集合Cp的长度，cq_j特征段集合Cq中第j个音频数据特征段，1≤j≤Lq，Lq为特征段集合Cq的长度，因为音频是无损音频，所以其每秒采样值数目比较大，通常都为每秒48000个采样点（即特征点），如音频p有8000个采样数据，音频q有6000个采样数据，则将p，q划分成长度为100个采样点的音频数据特征段，则音频p有80个音频数据特征段，每个音频数据段有100个特征点，音频q有60个音频数据特征段，每个音频数据特征段有100个特征点；

步骤2，利用CUDA（通用计算架构）方法，将音频p的特征段集合Cp和音频q的特征段集合Cq中每个特征段中的音频数据导入到GPU（Graphic Processing Unit，图形处理器，）上并行计算，得到音频p的特征段集合Cp的能量特征值序列Wp和音频q的特征段集合Cq的能量特征值序列Wq；

步骤3，将音频p的能量特征值序列Wp和音频q能量特征值序列Wq通过差值方法构建出一个Lp行Lq列能量矩阵Mat，其中Wp={wp₁，wp₂…wp_i…wp_Lp}，Wq={wq₁,wq₂…wq_j…wq_Lq}，即将能量特征值序列Wp中wp₁作为与Wq={wq₁,wq₂…wq_j…wq_Lq}中的每一个能量特征值做差值后作为一行，同理为wp₂，wp₃…wq_Lq中的每一个能量特征值都做差值得到对应行，组合后即得到Lp行Lq列能量矩阵Mat。

能量矩阵Mat可以采用如下所示矩阵结构：

$\mathrm{Mat}=\left[\begin{array}{cccccc}|w{p}_1-{\mathrm{wq}}_1|& |{\mathrm{wp}}_2-{\mathrm{wq}}_1|& \·\·\·& |{\mathrm{wp}}_i-{\mathrm{wq}}_1|& \·\·\·& |{\mathrm{wp}}_{\mathrm{Lp}}-{\mathrm{wq}}_1|\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ |{\mathrm{wp}}_1-{\mathrm{wq}}_j|& |{\mathrm{wp}}_2-{\mathrm{wq}}_j|& \·\·\·& |{\mathrm{wp}}_i-{\mathrm{wq}}_j|& \·\·\·& |{\mathrm{wp}}_{\mathrm{Lp}}-{\mathrm{wq}}_j|\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ |{\mathrm{wp}}_1-{\mathrm{wq}}_{\mathrm{Lq}}|& |{\mathrm{wp}}_2-{\mathrm{wq}}_{\mathrm{Lq}}|& \·\·\·& |{\mathrm{wp}}_i-{\mathrm{wq}}_{\mathrm{Lq}}|& \·\·\·& |{\mathrm{wp}}_{\mathrm{Lp}}-{\mathrm{wq}}_{\mathrm{Lq}}|\end{array}\right]$

；

步骤4，利用变形的后缀数组查找公共子串算法，查找出步骤3中所得的能量特征值序列Wp和能量特征值序列Wq中的公共部分，公共部分即是音频p和音频q相同区域段落的集合Seg；

步骤5，将能量矩阵Mat中元素数值小于设定阈值的区域连通，得到连通区域集合Vres；阈值设定为任意大于0小于1的实数，数值越小，音频比较的精度越高。

步骤6，将步骤5中所得的连通区域集合Vres作为结果，用波形表示；

步骤7，结束。

所述步骤2具体包括以下步骤：

步骤2a，将步骤1中所得到特征段cp_i中的特征点组成一个长度固定为N的向量D={P_i0,P_i1,P_i2…P_in…P_iN}，1≤i≤Lp，1≤n≤N，基于并将向量D通过CUDA（ComputeUnified DeviceArchitecture，通用计算架构）的cudamalloc函数加载到GPU内存中设置线程块block和线程thread，使得该设置动态最优的适合于向量D={P_i0,P_i1,P_i2…P_in…P_iN}；

步骤2b，将计算向量D={P_i0,P_i1,P_i2…P_in…P_iN}的特征值wp_i的过程分布到每个线程上，能量特征值wp_i的计算公式如下：

${\mathrm{wp}}_i=\frac{\underset{n=0}{\overset{N}{\mathrm{\Σ}}}[\mathrm{\α}\·\mathrm{energy}\left(P_{\mathrm{in}}\right)+\mathrm{\β}\·\mathrm{none}\left(P_{\mathrm{in}}\right)]}{\mathrm{dis}\left(N\right)}$

其中P_in为音频数据特征段cp_i上的特征点，N为音频数据特征段cp_i上的特征点的总数，0≤n≤N，α、β为设定的权值，α+β=1,0≤α≤1,0≤β≤1,α越大意味着该音频能量集中，β越大意味着该音频波动较大，energy(P_in)为音频数据特征段cp_i中特征点的能量信息，none(P_in)为音频数据特征段cp_i中特征点的过零点信息，dis(N)为音频数据特征段cp_i的距离信息；

步骤2c，将步骤1中所得到特征段cq_i中的特征点组成一个长度固定为N的向量D={P_j0,P_j1,P_j2…P_jn…P_jN}，1≤j≤Lq，1≤m≤N，基于并将向量D通过CUDA的cudamalloc函数加载到GPU内存中，合理设置线程块block和线程thread，使得该设置动态最优的适合于向量D={P_j0,P_j1,P_j2…P_jn…P_jN}；

步骤2d，将计算向量D={P_j0,P_j1,P_j2…P_jn…P_jN}的特征值wq_j的过程分布到每个线程上，特征值wq_i的计算公式如下：

${\mathrm{wq}}_j=\frac{\underset{m=0}{\overset{N}{\mathrm{\Σ}}}[\mathrm{\α}\·\mathrm{energy}\left(P_{\mathrm{jm}}\right)+\mathrm{\β}\·\mathrm{none}\left(P_{\mathrm{jm}}\right)]}{\mathrm{dis}\left(N\right)},$

其中P_jm为音频数据特征段cq_j上的特征点，N为音频数据特征段cq_j上的特征点的总数，0≤m≤N，由于音频数据特征段cp_i与音频数据特征段cq_j等长，因此两者的特征点的总数相等为N；α、β为设定的权值，α+β=1,0≤α≤1,0≤β≤1,α越大意味着该音频能量集中，β越大意味着该音频波动较大，energy(P_jm)为音频数据特征段cq_j中特征点的能量信息，none(P_jm)为音频数据特征段cq_j中特征点的过零点信息，dis(N)为音频数据特征段cq_j的距离信息。

步骤2e，步骤2b和步骤2d中得到的特征值数值范围不统一，将其归一到一个统一的尺度；

所述步骤4具体包括以下步骤：

步骤4a，将步骤3中的能量特征值序列Wp={wp₁，wp₂…wp_i…wp_Lp}和能量特征值序列Wq={wq₁,wq₂…wq_j…wq_Lq}(其中Lp、Lq分别为能量特征值序列的长度，1≤i≤Lp，1≤j≤Lq)结合起来，将能量特征值序列中数值设为一个字符，构建一个字符序列V_w={wp₁，wp₂…wp_i…wp_Lp,NA,wq₁,wq₂…wq_j…wq_Lq}，其中NA为隔断标记；

步骤4b，分别为步骤4a中Wp={wp₁，wp₂…wp_i…wp_Lp}和音频q能量特征值序列Wq={wq₁,wq₂…wq_j…wq_Lq}构建后缀数组，基本思路是计算Wp的所有后缀和Wq的所有后缀之间的最长公共前缀的长度，把最长公共前缀长度不小于k的部分全部加起来，k为设定的最小检索长度。扫描字符序列V_w{wp₁，wp₂…wp_i…wp_Lp,NA,wq₁,wq₂…wq_j…wq_Lq}，每遇到一个Wq的后缀就统计与前面的Wp的后缀能产生多少个长度不小于k的公共子串，这里Wp的后缀需要用一个单调的栈来高效的维护；

步骤4c，将步骤4b中返回的公共子串集合按照位置对应关系得出特征值向量Wp={wp₁，wp₂…wp_i…wp_Lp}和Wq={wq₁,wq₂…wq_j…wq_Lq}中对应的公共序列集合Seq；

步骤4d，将步骤2c特征值向量Wp={wp₁，wp₂…wp_i…wp_Lp}和Wq={wq₁,wq₂…wq_j…wq_Lq}中公共序列集合Seq合并整理，将公共交叉的数据段落合并，将连续的部分整理连接得到公共段落集合Seq^*；

所述步骤5具体包括以下步骤：

步骤5a，将步骤4中得到的相同区域的段落的集合Seg整理，对应于音频特征段得到音频特征段相同区域集合W_Seg={Seg₁(p,q),Seg₂(p,q)…Seg_h(p,q)…Seg_Lw(p,q)}，其中Seg_h(p,q)即为p，q两音频特征段的第h段公共区域，h介于1到Lw之间，Lw为音频特征段相同区域集合的长度；

步骤5b，按照特征段集合Cp={cp₁，cp₂…cp_i…cp_Lp}和Cp={cq₁，cq₂…cq_j…cq_Lq}将特征段相同区域集合W_Seg={Seg₁(p,q),Seg₂(p,q)…Seg_h(p,q)…Seg_Lw(p,q)}进行排序，遍历W_Seg，若Seg_s(p,q)和Seg_t(p,q)存在数据段上交叉，则将其合并，整理后得到新的特征段相同区域集合W_Seg ^*={Seg₁ ^*(p,q),Seg₂ ^*(p,q)…Seg_h ^*(p,q)…Seg_Lw ^*(p,q)}；

步骤5c，遍历步骤10b中W_Seg ^*={Seg₁ ^*(p,q),Seg₂ ^*(p,q)…Seg_h ^*(p,q)…Seg_Lw ^*(p,q)}，如满足下列条件则将其合并：

条件1：若存在Seg_s ^*(p,q)和Seg_t ^*(p,q)不相邻分别相隔1个特征段cp_x、cq_x，且Seg_s ^*(p,q)中音频p的特征段能量特征值为wp_s，音频q的能量特征段特征值为wq_s，Seg_t ^*(p,q)中音频p的特征段特征值为wp_t，音频q的特征段特征值为wq_t，特征段cp_x、cq_x特征值为wp_x、wq_x；

条件2：若wp_s=wq_s且wp_t=wq_t；

条件3；若wp_x=wq_x或|wp_x-wq_x|<T，T为阈值；

重复上述合并过程，直至不能合并为止，合并后得到的新的集合即为最简特征段相同区域集合W_Seg’={Seg₁’(p,q),Seg₂’(p,q)…Seg_h’(p,q)…Seg_Lw’(p,q)}；

步骤5d，验证步骤5c中得到最简特征段相同区域集合W_Seg’并对应到音频数据段得到连通区域集合Vres，则Vres即为音频中相同区域的集合；

所述步骤2b具体包括以下步骤：

步骤2ba，计算步骤2a中的向量D={P_i0,P_i1,P_i2…P_in…P_iN}中特征点P_in的能量信息，记为energy，则P_in的能量信息为energy(P_in)，energy(P_in)的计算公式如下：

$e\left(P_{\mathrm{in}}\right)={|P_{\mathrm{in}}-\mathrm{Avg}\left(\underset{n=0}{\overset{N}{\mathrm{\&Sigma;}}}{p}_{\mathrm{in}}\right)|}^2$

其中P_in为向量D={P_i0,P_i1,P_i2…P_in…P_iN}第n点的数据，1≤i≤Lp，1≤n≤N，

为该特征段的平均数据值；

步骤2bb，计算步骤2a中的向量D={P_i0,P_i1,P_i2…P_in…P_iN}中特征点P_in的过零点信息，记为none，则P_in的过零点信息为none(P_in)，none(P_in)的计算公式如下：

$\mathrm{none}\left(P_{\mathrm{in}}\right)=\left\{\begin{array}{c}0,P_{\mathrm{in}}\&NotEqual;0\\ 1,P_{\mathrm{in}}=0\end{array}\right.$

其中P_in为向量D={P_i0,P_i1,P_i2…P_in…P_iN}第n点的数据，1≤i≤Lp，1≤n≤N；

步骤2bc，根据音频的特征性质动态地设置权值α、β，使得能量信息和过零点信息按照比重影响特征值；

步骤2bd，计算步骤2a中的向量D={P_i0,P_i1,P_i2…P_in…P_iN}的能量距离dis(N)，dis(N)的计算公式如下：

dis(N)=γN

其中γ为与能量信息相关的参数，N为向量D={P_i0,P_i1,P_i2…P_in…P_iN}的长度；

步骤2be，根据步骤2ba-2bd中每个P_in的能量信息energy(P_in)，过零点信息none(P_in)，能量距离dis(N)以及权值α、β按照步骤2b中的公式计算出特征段cp_i的特征值wp_i；

所述步骤2d具体包括以下步骤：

步骤2da，计算步骤2a中的向量D={P_j0,P_j1,P_j2…P_jm…P_jN}中特征点P_jm的能量信息，记为energy，则P_jm的能量信息为energy(P_jm)，energy(P_jm)的计算公式如下：

$e\left(P_{\mathrm{jm}}\right)={|P_{\mathrm{jm}}-\mathrm{Avg}\left(\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}{p}_{\mathrm{jm}}\right)|}^2$

其中P_jm为向量D={P_j0,P_j1,P_j2…P_jm…P_jN}第m点的数据，1≤i≤Lp，1≤m≤N，

为该特征段的平均数据值；

步骤2db，计算步骤2a中的向量D={P_j0,P_j1,P_j2…P_jm…P_jN}中特征点P_jm的过零点信息，记为none，则P_jm的过零点信息为none(P_jm)，none(P_jm)的计算公式如下：

$\mathrm{none}\left(P_{\mathrm{jm}}\right)=\left\{\begin{array}{c}0,P_{\mathrm{jm}}\&NotEqual;0\\ 1,P_{\mathrm{jm}}=0\end{array}\right.$

其中P_jm为向量D={P_j0,P_j1,P_j2…P_jm…P_jN}第m点的数据，1≤i≤Lp，1≤m≤N；

步骤2dc，根据音频的特征性质动态地设置权值α、β，使得能量信息和过零点信息按照比重影响特征值；

步骤2dd，计算步骤2a中的向量D={P_j0,P_j1,P_j2…P_jm…P_jN}的能量距离dis(N)，dis(N)的计算公式如下：

dis(N)=γN

其中γ为与能量信息相关的参数，0≤γ≤1，N为向量D={P_j0,P_j1,P_j2…P_jm…P_jN}的长度；

步骤2de，根据步骤2ba-2bd中每个P_in的能量信息energy(P_in)，过零点信息none(P_in)，能量距离dis(N)以及权值α、β按照步骤2d中的公式计算出特征段cq_j的特征值wq_j；

所述步骤4b具体包括以下步骤：

步骤4ba，为步骤4a中的特征值向量Wp={wp₁，wp₂…wp_i…wp_Lp}建立后缀数组S_Wp，它保存1-Lp的某个排列，S_Wp={S_Wp1,S_Wp2…S_Wpi…S_WpLp}，保证suffix(wp_[i])<suffix(wp_[i+1])suffix为从i位置开始到整个序列结尾的一个子序列，其中suffix(wp_[i])=Wp_[i-Lp]；同样为特征值向量Wq={wq₁，wq₂…wq_j…wq_Lp}建立其后缀数组S_Wq={S_Wq1,S_Wq2…S_Wqj…S_WqLq}；

步骤4bb，计算步骤4ba中的后缀数组S_Wp={S_Wp1,S_Wp2…S_Wpi…S_WpLp}计算其对应位置的height值，即得到一个名次数组Rank_Wp，Rank_Wpi的计算公式如下：

Rank_Wpi=No(suffix(wp_i))

其中No为排名函数，suffix(wp_i)为对应i位置的后缀数组；同样为步骤4ba中的后缀数组S_Wq={S_Wq1,S_Wq2…S_Wqj…S_WqLq}计算名次数组Rank_Wq；

步骤4bc，按照Rank_Wp和Rank_Wq的值分组，快速统计每组中后缀之间的最长公共前缀之和，求两个后缀的最长公共前缀可以转化为某个区间上的最小值问题，即RMQ(Range Minimum Query)问题；

步骤4bd，扫描步骤4a中的V_w{wp₁,wp₂…wp_i…wp_Lp,NA,wq₁,wq₂…wq_j…wq_Lq}，每遇到一个S_Wqj后缀就统计与前面的S_Wp={S_Wp1,S_Wp2…S_Wpi…S_WpLp}的后缀能产生多少个长度不小于k的公共子串，k为设定的最小检索长度，这里S_Wp={S_Wp1,S_Wp2…S_Wpi…S_WpLp}的后缀需要用一个单调的栈来高效的维护；

步骤4be，计算S_Wp={S_Wp1,S_Wp2…S_Wpi…S_WpLp}的所有suffix(wp_i)后缀和S_Wq={S_Wq1,S_Wq2…S_Wqj…S_WqLq}的所有suffix(wq_j)后缀之间的最长公共前缀的长度，把最长公共前缀长度不小于k的部分全部加起来，得到一个公共子序列集合；

有益效果：本发明的显著优点是比较方法具有很好的效率并且比较结果具有很高的精度。具体而言，本发明将分别从算法的精度和效率方面来阐述：（1）算法的精度保证：采用多重特征来对音频文件提取特征，从音频采样数据，音频波形性质以及能量矩阵的方法多尺度的提取特征，提高比较精度。（2）算法的效率保证：本发明采用了GPU上的CUDA大大加速了比较速度，将原本的NP问题在短时间实时解决，可以很好满足工业级别的运用。因此，本发明能更好地保证检索的效率，满足实时响应的需求。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述或其他方面的优点将会变得更加清楚。

图1为本发明方法的流程图。

图2为本发明实施例的流程图。

图3变形的后缀数组查找公共子串算法的流程图。

具体实施方式：

如图1所示，音频显示和检索装置从音频存储设备获取数字音频，假设音频存储设备存储了若干段音频，装置同时接受Num(Num≤16)段音频，然后装置生成音频的波形特征显示，并提供播放功能，用户在Num段音频中任选两段音频作为比较对象，装置生成音频的特征表示，使用采用多重特征(音频采样数据，音频波形性质以及能量矩阵)表示，这样每段音频有一个特征值向量表示，通过后缀数组算法得出比较结果并验证。

如图2所示，音频快速比较部分包括以下步骤：

步骤1，将待比较的音频p和音频q，分别按照其音频的长度划分出长度相等的两组音频数据特征段，其中音频p划分为特征段集合Cp，Cp={cp₁，cp₂…cp_i…cp_Lp}，音频q划分为特征段集合Cq，Cq={cq₁，cq₂…cq_j…cq_Lq}，其中cp_i表示特征段集合Cp中第i个音频数据特征段，1≤i≤Lp，Lp为特征段集合Cp的长度，cq_j特征段集合Cq中第j个音频数据特征段，1≤j≤Lq，Lq为特征段集合Cq的长度，因为音频是无损音频，所以其每秒采样值数目比较大，通常都为每秒48000个采样点（即特征点），如音频p有8000个采样数据，音频q有6000个采样数据，则将p，q划分成长度为100个采样点的音频数据特征段，则音频p有80个音频数据特征段，每个音频数据段有100个特征点，音频q有60个音频数据特征段，每个音频数据特征段有100个特征点；

步骤2，将音频p的特征段集合Cp和音频q的特征段集合Cq中每个特征段中的音频数据导入到GPU上并行计算，得到音频p的特征段集合Cp的能量特征值序列Wp和音频q的特征段集合Cq的能量特征值序列Wq；

步骤3，将音频p的能量特征值序列Wp和音频q能量特征值序列Wq通过差值方法构建出一个Lp行Lq列能量矩阵Mat，其中Wp={wp₁，wp₂…wp_i…wp_Lp}，Wq={wq₁,wq₂…wq_j…wq_Lq}，即将Wp中wp₁作为与Wq={wq₁,wq₂…wq_j…wq_Lq}中的每一个能量特征值做差值后作为一行，同理为wp₂，wp₃…wq_Lq中的每一个能量特征值都做差值得到对应行，组合后即得到Lp行Lq列能量矩阵Mat，Mat如下所示：

$\mathrm{Mat}=\left[\begin{array}{cccccc}|w{p}_1-{\mathrm{wq}}_1|& |{\mathrm{wp}}_2-{\mathrm{wq}}_1|& \&CenterDot;\&CenterDot;\&CenterDot;& |{\mathrm{wp}}_i-{\mathrm{wq}}_1|& \&CenterDot;\&CenterDot;\&CenterDot;& |{\mathrm{wp}}_{\mathrm{Lp}}-{\mathrm{wq}}_1|\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ |{\mathrm{wp}}_1-{\mathrm{wq}}_j|& |{\mathrm{wp}}_2-{\mathrm{wq}}_j|& \&CenterDot;\&CenterDot;\&CenterDot;& |{\mathrm{wp}}_i-{\mathrm{wq}}_j|& \&CenterDot;\&CenterDot;\&CenterDot;& |{\mathrm{wp}}_{\mathrm{Lp}}-{\mathrm{wq}}_j|\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ |{\mathrm{wp}}_1-{\mathrm{wq}}_{\mathrm{Lq}}|& |{\mathrm{wp}}_2-{\mathrm{wq}}_{\mathrm{Lq}}|& \&CenterDot;\&CenterDot;\&CenterDot;& |{\mathrm{wp}}_i-{\mathrm{wq}}_{\mathrm{Lq}}|& \&CenterDot;\&CenterDot;\&CenterDot;& |{\mathrm{wp}}_{\mathrm{Lp}}-{\mathrm{wq}}_{\mathrm{Lq}}|\end{array}\right]$

步骤4，利用变形的后缀数组查找公共子串算法，查找出步骤3中所得的能量特征值序列Wp和能量特征值序列Wq中的公共部分，公共部分即是音频p和音频q相同区域段落的集合Seg；

步骤5，将能量矩阵Mat中元素数值小于设定阈值的区域连通，得到连通区域集合Vres，阈值取值范围为0~1的实数，数值越小，音频比较的精度越高；

步骤6，将步骤5中所得的连通区域集合Vres作为结果，用波形表示；

步骤7，结束。

所述步骤2具体包括以下步骤：

步骤2a，将步骤1中所得到特征段cp_i中的特征点组成一个长度固定为N的向量D={P_i0,P_i1,P_i2…P_in…P_iN}，1≤i≤Lp，1≤n≤N，基于并将向量D通过CUDA的cudamalloc函数加载到GPU内存中，合理设置线程块block和线程thread，使得该设置动态最优的适合于向量D={P_i0,P_i1,P_i2…P_in…P_iN}；

步骤2b，将计算向量D={P_i0,P_i1,P_i2…P_in…P_iN}的特征值wp_i的过程分布到每个线程上，特征值wp_i的计算公式如下：

${\mathrm{wp}}_i=\frac{\underset{n=0}{\overset{N}{\mathrm{\&Sigma;}}}[\mathrm{\&alpha;}\&CenterDot;\mathrm{energy}\left(P_{\mathrm{in}}\right)+\mathrm{\&beta;}\&CenterDot;\mathrm{none}\left(P_{\mathrm{in}}\right)]}{\mathrm{dis}\left(N\right)}$

其中P_in为音频数据特征段cp_i上的特征点，N为音频数据特征段cp_i上的特征点的总数，0≤n≤N，α、β为设定的权值，α+β=1,0≤α≤1,0≤β≤1,α越大意味着该音频能量集中，β越大意味着该音频波动较大，energy(P_in)为音频数据特征段cp_i中特征点的能量信息，none(P_in)为音频数据特征段cp_i中特征点的过零点信息，dis(N)为音频数据特征段cp_i的距离信息；

步骤2c，将步骤1中所得到特征段cq_i中的特征点组成一个长度固定为N的向量D={P_j0,P_j1,P_j2…P_jn…P_jN}，1≤j≤Lq，1≤m≤N，基于并将向量D通过CUDA的cudamalloc函数加载到GPU内存中，合理设置线程块block和线程thread，使得该设置动态最优的适合于向量D={P_j0,P_j1,P_j2…P_jn…P_jN}；

步骤2d，将计算向量D={P_j0,P_j1,P_j2…P_jn…P_jN}的特征值wq_j的过程分布到每个线程上，特征值wq_i的计算公式如下：

${\mathrm{wq}}_j=\frac{\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}[\mathrm{\&alpha;}\&CenterDot;\mathrm{energy}\left(P_{\mathrm{jm}}\right)+\mathrm{\&beta;}\&CenterDot;\mathrm{none}\left(P_{\mathrm{jm}}\right)]}{\mathrm{dis}\left(N\right)},$

其中P_jm为音频数据特征段cq_j上的特征点，N为音频数据特征段cq_j上的特征点的总数，0≤m≤N，α、β为设定的权值，α+β=1,0≤α≤1,0≤β≤1,α越大意味着该音频能量集中，β越大意味着该音频波动较大，energy(P_jm)为音频数据特征段cq_j中特征点的能量信息，none(P_jm)为音频数据特征段cq_j中特征点的过零点信息，dis(N)为音频数据特征段cq_j的距离信息。

步骤2e，步骤2b和步骤2d中得到的特征值数值范围不统一，我们将其归一到一个统一的尺度；

所述步骤4具体包括以下步骤：

步骤4a，将步骤3中的能量特征值序列Wp={wp₁，wp₂…wp_i…wp_Lp}和能量特征值序列Wq={wq₁,wq₂…wq_j…wq_Lq}(其中Lp、Lq分别为能量特征值序列的长度，1≤i≤Lp，1≤j≤Lq)结合起来，将能量特征值序列中数值设为一个字符，构建一个字符序列V_w={wp₁，wp₂…wp_i…wp_Lp,NA,wq₁,wq₂…wq_j…wq_Lq}，其中NA为隔断标记；

步骤4b，分别为步骤4a中能量特征值序列Wp={wp₁，wp₂…wp_i…wp_Lp}和音频q能量特征值序列Wq={wq₁,wq₂…wq_j…wq_Lq}构建后缀数组，基本思路是计算Wp的所有后缀和Wq的所有后缀之间的最长公共前缀的长度，把最长公共前缀长度不小于k的部分全部加起来，k为设定的最小检索长度。扫描字符序列V_w={wp₁，wp₂…wp_i…wp_Lp,NA,wq₁,wq₂…wq_j…wq_Lq}，每遇到一个Wq的后缀就统计与前面的Wp的后缀能产生多少个长度不小于k的公共子串，这里Wp的后缀需要用一个单调的栈来高效的维护；

步骤4c，将步骤4b中返回的公共子串集合按照位置对应关系得出特征值向量Wp={wp₁，wp₂…wp_i…wp_Lp}和Wq={wq₁,wq₂…wq_j…wq_Lq}中对应的公共序列集合Seq；

步骤4d，将步骤2c特征值向量Wp={wp₁，wp₂…wp_i…wp_Lp}和Wq={wq₁,wq₂…wq_j…wq_Lq}中公共序列集合Seq合并整理，将公共交叉的数据段落合并，将连续的部分整理连接得到公共段落集合Seq^*；

所述步骤5具体包括以下步骤：

步骤5a，将步骤4中得到的相同区域的段落的集合Seg整理，对应于音频特征段得到音频特征段相同区域集合W_Seg={Seg₁(p,q),Seg₂(p,q)…Seg_h(p,q)…Seg_Lw(p,q)}，其中Seg_h(p,q)即为p，q两音频特征段的第h段公共区域，h介于1到Lw之间，Lw为音频特征段相同区域集合的长度；

步骤5b，按照特征段集合Cp={cp₁，cp₂…cp_i…cp_Lp}和Cp={cq₁，cq₂…cq_j…cq_Lq}将特征段相同区域集合W_Seg={Seg₁(p,q),Seg₂(p,q)…Seg_h(p,q)…Seg_Lw(p,q)}进行排序，遍历W_Seg，若Seg_s(p,q)和Seg_t(p,q)存在数据段上交叉，则将其合并，整理后得到新的特征段相同区域集合W_Seg ^*={Seg₁ ^*(p,q),Seg₂ ^*(p,q)…Seg_h ^*(p,q)…Seg_Lw ^*(p,q)}；

步骤5c，遍历步骤10b中W_Seg ^*={Seg₁ ^*(p,q),Seg₂ ^*(p,q)…Seg_h ^*(p,q)…Seg_Lw ^*(p,q)}，如满足下列条件则将其合并：

条件1：若存在Seg_s ^*(p,q)和Seg_t ^*(p,q)不相邻分别相隔1个特征段cp_x、cq_x，且Seg_s ^*(p,q)中音频p的特征段能量特征值为wp_s，音频q的能量特征段特征值为wq_s，Seg_t ^*(p,q)中音频p的特征段特征值为wp_t，音频q的特征段特征值为wq_t，特征段cp_x、cq_x特征值为wp_x、wq_x；

条件2：若wp_s=wq_s且wp_t=wq_t；

条件3；若wp_x=wq_x或|wp_x-wq_x|<T，T为阈值；

重复上述合并过程，直至不能合并为止，合并后得到的新的集合即为最简特征段相同区域集合W_Seg’={Seg₁’(p,q),Seg₂’(p,q)…Seg_h’(p,q)…Seg_Lw’(p,q)}；

步骤5d，验证步骤5c中得到最简特征段相同区域集合W_Seg’并对应到音频数据段得到连通区域集合Vres，则Vres即为音频中相同区域的集合；

所述步骤2b具体包括以下步骤：

步骤2ba，计算步骤2a中的向量D={P_i0,P_i1,P_i2…P_in…P_iN}中特征点P_in的能量信息，记为energy，则P_in的能量信息为energy(P_in)，energy(P_in)的计算公式如下：

$e\left(P_{\mathrm{in}}\right)={|P_{\mathrm{in}}-\mathrm{Avg}\left(\underset{n=0}{\overset{N}{\mathrm{\&Sigma;}}}{p}_{\mathrm{in}}\right)|}^2$

其中P_in为向量D={P_i0,P_i1,P_i2…P_in…P_iN}第n点的数据，1≤i≤Lp，1≤n≤N，为该特征段的平均数据值；

步骤2bb，计算步骤2a中的向量D={P_i0,P_i1,P_i2…P_in…P_iN}中特征点P_in的过零点信息，记为none，则P_in的过零点信息为none(P_in)，none(P_in)的计算公式如下：

$\mathrm{none}\left(P_{\mathrm{in}}\right)=\left\{\begin{array}{c}0,P_{\mathrm{in}}\&NotEqual;0\\ 1,P_{\mathrm{in}}=0\end{array}\right.$

其中P_in为向量D={P_i0,P_i1,P_i2…P_in…P_iN}第n点的数据，1≤i≤Lp，1≤n≤N；

步骤2bc，根据音频的特征性质动态地设置权值α、β，使得能量信息和过零点信息按照比重影响特征值；

步骤2bd，计算步骤2a中的向量D={P_i0,P_i1,P_i2…P_in…P_iN}的能量距离dis(N)，dis(N)的计算公式如下：

dis(N)=γN

其中γ为与能量信息相关的参数，N为向量D={P_i0,P_i1,P_i2…P_in…P_iN}的长度；

步骤2be，根据步骤2ba-2bd中每个P_in的能量信息energy(P_in)，过零点信息none(P_in)，能量距离dis(N)以及权值α、β按照步骤2b中的公式计算出特征段cp_i的特征值wp_i；

所述步骤2d具体包括以下步骤：

步骤2da，计算步骤2a中的向量D={P_j0,P_j1,P_j2…P_jm…P_jN}中特征点P_jm的能量信息，记为energy，则P_jm的能量信息为energy(P_jm)，energy(P_jm)的计算公式如下：

$e\left(P_{\mathrm{jm}}\right)={|P_{\mathrm{jm}}-\mathrm{Avg}\left(\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}{p}_{\mathrm{jm}}\right)|}^2$

其中P_jm为向量D={P_j0,P_j1,P_j2…P_jm…P_jN}第m点的数据，1≤i≤Lp，1≤m≤N，

为该特征段的平均数据值；

步骤2db，计算步骤2a中的向量D={P_j0,P_j1,P_j2…P_jm…P_jN}中特征点P_jm的过零点信息，记为none，则P_jm的过零点信息为none(P_jm)，none(P_jm)的计算公式如下：

$\mathrm{none}\left(P_{\mathrm{jm}}\right)=\left\{\begin{array}{c}0,P_{\mathrm{jm}}\&NotEqual;0\\ 1,P_{\mathrm{jm}}=0\end{array}\right.$

其中P_jm为向量D={P_j0,P_j1,P_j2…P_jm…P_jN}第m点的数据，1≤i≤Lp，1≤m≤N；

步骤2dc，根据音频的特征性质动态地设置权值α、β，使得能量信息和过零点信息按照比重影响特征值；

步骤2dd，计算步骤2a中的向量D={P_j0,P_j1,P_j2…P_jm…P_jN}的能量距离dis(N)，dis(N)的计算公式如下：

dis(N)=γN

其中γ为与能量信息相关的参数，N为向量D={P_j0,P_j1,P_j2…P_jm…P_jN}的长度；

步骤2de，根据步骤2ba-2bd中每个P_in的能量信息energy(P_in)，过零点信息none(P_in)，能量距离dis(N)以及权值α、β按照步骤2d中的公式计算出特征段cq_j的特征值wq_j；

所述步骤4b具体包括以下步骤：

步骤4ba，为步骤4a中的特征值向量Wp={wp₁，wp₂…wp_i…wp_Lp}建立后缀数组S_Wp，它保存1-Lp的某个排列，S_Wp={S_Wp1,S_Wp2…S_Wpi…S_WpLp}，保证suffix(wp_[i])<suffix(wp_[i+1])suffix为从i位置开始到整个序列结尾的一个子序列，其中suffix(wp_[i])=Wp_[i-Lp]；同样为特征值向量Wq={wq₁，wq₂…wq_j…wq_Lp}建立其后缀数组S_Wq={S_Wq1,S_Wq2…S_Wqj…S_WqLq}；

步骤4bb，计算步骤4ba中的后缀数组S_Wp={S_Wp1,S_Wp2…S_Wpi…S_WpLp}计算其对应位置的height值，即得到一个名次数组Rank_Wp，Rank_Wpi的计算公式如下：

Rank_Wpi=No(suffix(wp_i))

其中No为排名函数，suffix(wp_i)为对应i位置的后缀数组；同样为步骤4ba中的后缀数组S_Wq={S_Wq1,S_Wq2…S_Wqj…S_WqLq}计算名次数组Rank_Wq；

步骤4bc，按照Rank_Wp和Rank_Wq的值分组，快速统计每组中后缀之间的最长公共前缀之和，求两个后缀的最长公共前缀可以转化为某个区间上的最小值问题，即RMQ(Range Minimum Query)问题；

步骤4bd，扫描步骤4a中的V_w{wp₁,wp₂…wp_i…wp_Lp,NA,wq₁,wq₂…wq_j…wq_Lq}，每遇到一个S_Wqj后缀就统计与前面的S_Wp={S_Wp1,S_Wp2…S_Wpi…S_WpLp}的后缀能产生多少个长度不小于k的公共子串，k为设定的最小检索长度，这里S_Wp={S_Wp1,S_Wp2…S_Wpi…S_WpLp}的后缀需要用一个单调的栈来高效的维护；

步骤4be，计算S_Wp={S_Wp1,S_Wp2…S_Wpi…S_WpLp}的所有suffix(wp_i)后缀和S_Wq={S_Wq1,S_Wq2…S_Wqj…S_WqLq}的所有suffix(wq_j)后缀之间的最长公共前缀的长度，把最长公共前缀长度不小于k的部分全部加起来，得到一个公共子序列集合；如图3所示，给出了用本发明的方法进行相关测试的例子。用于测试的硬件环境是：Intel-Xeon 56202.66GHz16G内存和Nvidia GTX440显卡。软件环境是QT Creator、Centos4.5以及CUDAtoolkit4.0。本发明用C++语言实现了本发明提出的方法。测试音频来源于Disney。本发明从中抽取了1组音频数据进行了测试。由图3可见，本发明的方法具有很高的比较精度。

本发明提供了一种音频比较方法的思路及方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (2)

1.一种音频比较方法，其特征在于包括以下步骤：

步骤1，将待比较的音频p和音频q，分别按照其音频的长度划分出长度相等的两组音频数据特征段，其中音频p划分为特征段集合Cp，Cp={cp₁，cp₂…cp_i…cp_Lp}，音频q划分为特征段集合Cq，Cq={cq₁，cq₂…cq_j…cq_Lq}，其中cp_i表示特征段集合Cp中第i个音频数据特征段，1≤i≤Lp，Lp为特征段集合Cp的长度，cq_j特征段集合Cq中第j个音频数据特征段，1≤j≤Lq，Lq为特征段集合Cq的长度；

步骤2，将音频p的特征段集合Cp和音频q的特征段集合Cq中每个特征段中的音频数据导入到GPU上并行计算，得到音频p的特征段集合Cp的能量特征值序列Wp和音频q的特征段集合Cq的能量特征值序列Wq；

步骤3，将音频p的能量特征值序列Wp和音频q能量特征值序列Wq通过差值方法构建出一个Lp行Lq列能量矩阵Mat；

步骤4，利用变形的后缀数组查找公共子串算法，查找出步骤3中所得的能量特征值序列Wp和能量特征值序列Wq中的公共部分，公共部分即是音频p和音频q相同区域段落的集合Seg；

步骤5，将能量矩阵Mat中元素数值小于设定阈值的区域连通，得到连通区域集合Vres，阈值取值范围为0～1的实数；

步骤6，将步骤5中所得的连通区域集合Vres作为结果，用波形表示；

步骤7，结束；

在GPU上使用通用计算架构CUDA方法计算音频数据特征段cp_i的能量特征值wp_i和音频数据特征段cq_j的能量特征值wq_j，并将能量特征值wpi和能量特征值wq_j分布到GPU的线程上：

能量特征值wp_i的计算方法如下：

$w{p}_i=\frac{\underset{n=0}{\overset{N}{\mathrm{\&Sigma;}}}[\mathrm{\&alpha;}\&CenterDot;\mathrm{energy}\left(P_{\mathrm{in}}\right)+\mathrm{\&beta;}\&CenterDot;\mathrm{none}\left(P_{\mathrm{in}}\right)]}{\mathrm{dis}\left(N\right)},$

其中P_in为音频数据特征段cp_i上的特征点，N为音频数据特征段cp_i上的特征点的总数，0≤n≤N，α、β为设定的权值，α+β=1,0≤α≤1,0≤β≤1,energy(P_in)为音频数据特征段cp _i中特征点的能量信息，none(P_in)为音频数据特征段cp_i中特征点的过零点信息，dis(N)为音频数据特征段cp_i的距离信息；

能量特征值wq_j的计算方法如下：

$w{p}_j=\frac{\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}[\mathrm{\&alpha;}\&CenterDot;\mathrm{energy}\left(P_{\mathrm{jm}}\right)+\mathrm{\&beta;}\&CenterDot;\mathrm{none}\left(P_{\mathrm{jm}}\right)]}{\mathrm{dis}\left(N\right)},$

其中P_jm为音频数据特征段cq_j上的特征点，N为音频数据特征段cq_j上的特征点的总数，0≤m≤N，α、β为设定的权值，α+β=1,0≤α≤1,0≤β≤1,energy(P_jm)为音频数据特征段cq_j中特征点的能量信息，none(P_jm)为音频数据特征段cq_j中特征点的过零点信息，dis(N)为音频数据特征段cq_j的距离信息；

步骤5中，相同连通区域的判定标准为：音频数据特征段的能量特征值同时满足以下两个条件，则音频p的能量特征值序列Wp和音频q能量特征值序列Wq中第s、s+1、s+2三段数据是相同的并且可连通：

条件1：wp_(s)=wq_(s)且wp_(s+2)=wq_(s+2)；

条件2：wp_(s+1)=wq_(s+1)或者|wp_(s+1)－wq_(s+1)|≤T，T为阈值；

其中wp_(s)是音频p的能量特征值序列Wp的第s、s+1、s+2段数据，其中wq_(s)是音频q的能量特征值序列Wq的第s、s+1、s+2段数据。

2.根据权利要求1所述的一种音频比较方法，其特征在于，步骤4中，利用变形的后缀数组查找公共子串的算法，算法中查找的时间复杂度为O(nlgn)，将音频数据特征段的能量特征值作为变形的后缀数组查找公共子串的算法的字符串中的一个字符，得到若干段长度大于k的字符串，k为设定的最小区别长度，将字符串对应的位置转换为音频数据特征段的位置，得出比较结果。

Images