声纹识别

声纹识别系统

摘要

本文首先通过用层次分析法（AHP ）构建了影响声纹识别的八大因素，并将其进行量化处理，得到了合理的权重比。同时构建Mel 倒谱系统提取说话人的主要特征倒频谱（MFCC 参数），利用提取的 MFCC 参数训练话者的 GMM （高斯混合模型），得到专属某话者的 GMM 声纹模型。接着通过语音数据对构建的声纹识别模型进行评价。发现得到的模型虽然可以接受，但其准确率不是很高。为了使模型的精确率提高，我们通过利用改进的K-means 算法来将模型进行调整。

问题一：

我们通过层次分析发得出语者的声学特征在语音识别中所占的权重最大，而且限于目前的技术，我们最终将能描绘语者声道结构的Mel 倒频谱系数作为语音的特征向量，该方法能很好的模拟人耳对不同频率的感知特性，具有很好的稳定性和准确度。然后我们利用matlab 编程将提供的语音数据的MFCC 系数提取出来作为声纹识别系统的训练样本。

问题二：基于高斯混合模型(GMM),我们用EM 算法来估计GMM 中的未知参数，由K-means 算法来得到EM 算法的初始值，通过不断地迭代EM 算法更新GMM 模型中的权重（w ），均值(u)和方差矩阵,直到最终使

()*()P X P X λλε-<，ε取10-5。

问题三：

根据问题一的MFCC 模型和问题二的GMM 模型，利用matlab 编程，用样本语音数据去训练GMM 模型得到样本库，让后用未知的测试语音样本利用模式匹配法去与样本语音库进行匹配，得到一系列的匹配概率，概率最大者对应的语者为测试语音的主人。我们得出的结论是EM 算法的初值对最终的识别率有很大影响，用来确定GMM 初始值的K-means 算法有比较大的误差，需要改进。

问题四：

基于问题三的测试和查阅的相关资料，我们认为传统的K-means 算法忽略特征矢量各维分量对识别的影响，因而无法得到令人满意的识别结果。由数理统计知，方差是用来衡量样本数据的离散程度的。如果特征矢量的分布稀疏，则它的方差就会很大，反之，那些分布较集中的矢量方差比较小。为此，本文提出了基于方差的加权几何距离，在聚类过程中，对特征矢量的各分量按方差大小进行加权，加权因子为矢量方差的倒数。

关键字：声纹识别 层次分析法（AHP ） 美尓倒频谱系数（MFCC ） 高斯混合模型（GMM ） 期望值最大化算法（EM ） K-means 算法

一、问题重述

生物认证是通过人体特征进行人物身份认证的重要手段，包括人脸识别、虹膜识别、指纹识别、体态识别和声纹识别等方法，它们被广泛地应用于国防、军事和民用等诸多领域，推动着信息产业的蓬勃发展。

声纹识别又称说话人识别，是通过分析说话人的语音的特点，利用计算机自动识别说话人的身份。声纹识别又分为文本相关和文本无关说话人的识别两大类。文本无关是指说话人的说话内容不受限制，文本相关是指说话人的说话内容必须是指定的内容。

评价声纹识别方法或模型的最重要指标是识别准确率，准确率越高越宜于应用于实际。评价的另外一个重要指标就是识别所花的时间，即输入待识别的说话人语音到识别结果的输出说话的时间，该时间越短约好。

请试着利用采集到的说话人语音数据（文本相关和文本无关两类），完成以下几点：

1）建立代表说话人身份的语音特征模型，通过编程从语音数据中提取说话人特征；

2）利用说话人特征建立说话人模型或声纹识别模型；

3）利用语音数据评价你的说话人模型或声纹识别模型评价的好坏；

4）分析影响准确率的因素，采取措施进一步提高声纹识别的准确率。

二、问题分析

问题一：

在没有人的主观倾向影响的情况下，人类的听觉系统其实就是比较好的说话人识别系统，具有很高的准确性和抗干扰性。虽然人的听觉系统分辨说话人的机理不一定是最佳的声纹识别方法，但是在目前的技术下，如果能达到人类说话人的水平，也是相当可观的。MFCC是基于短时傅里叶变换谱分析，并考虑人类的听觉特性后得到的一种倒频特征。心理学的研究进一步表明，人类对于声音音调的感觉都不是线性的，这必然导致人们去定义新的频率单位。符合人的听觉特性的频率则以Mel或Bark作为单位；

问题二：就声学观点而言，每一位语者声学特征能分出很多类别，诸如：鼻音，摩擦音，元音等，类别还可能更细致，而且每一个声学类别下又有广阔的语音空间，查阅资料发现，高斯混合模型(GMM)能好的二者兼顾，于是可以建立GMM模型对每位语者的声学特征进行仿真。

问题三：基于上面两问的模型，用matlab编程对实际的数据进行识别率的分析即可；

问题四：

传统的K-means算法在聚类过程中采取几何距离就近原则，这只有当类的自然布为球状或接近于球状时，也就是说当类中各矢量的分布接近于相等时才能取得识别效果。而语音信号是近似椭圆状的正态分布，各矢量的分布不等，若用传统K．means算法作为GMM 的初始化方法，将会忽略特征矢量各维分量对识别的影响，无法得到令人满意的识别结果。

为此，对传统K。means算法进行了改进，语音信号各矢量的分布不同，且对于不同的说话人这种差别的程度又不一样。好的特征矢量应该对同类语音的离散程度小，对不同类语音的离散程度大。由数理统计知，方差是用来衡量样本数据的离散程度的。如果特征矢量的分布稀疏，则它的方差就会很大，反之，那些分布较集中的矢量方差比较小。为此，本文提出了基于方差的加权几何距离，在聚类过程中，对特征矢量的各分量按方差大小进行加权，加权因子为矢量方差的倒数。

三、模型假设与符号说明

3.1 模型假设：

（1）本文选用文本无关的声纹识别文件；

（2）假设提供的数据完整、可靠、有效；

（3）假设查阅的资料和一些结论正确可靠；

3.2 符号说明：

()

x n：预处理后的语音信号。

()

X k：将预处理后的语音信号进行短时傅里叶变换后得到的频谱。

Q：频率滤波器个数。

()

S q：作对数转换后的Mel频谱。

()

C n：Mel频谱倒谱系数。

四、模型的建立及求解

4.1 问题一模型的建立与求解：

4.1.1 层次结构模型建立

图1 层次结构示意图

首先分析准则层对目标层的影响，设有n 个因素（x 1,x 2,…x n ）,用a ij 表示x i 和x j 对上层目标的影响比。

相对重要程度：

a ij

定义

1 若 xi 等价于xj: 赋值 1 3 若 xi 比xj 重要 : 赋值 3 5 若 xi 比xj 重要得多: 赋值 5 7 若 xi 远远比xj 重要: 赋值 7

9 若xi 是最重要的

2,4,6,8 重要程度等级介于 xi 和 xj 之间 1/2,1/3,..,1/9

对应于以上等级的xi 和xj 之间的关系

表1 相对重要程度a ij 取值情况

提取说话人特

征

沙哑音 鼻音 语音习惯

个人特点 声学特征 倒频谱

频谱

基音

共振峰 词法特征

韵律特征

口音

通道信息

由网络资源查询可对比声学特征，语音习惯，个人特点，鼻音，沙哑音五个元素得到以下系数矩阵：

ij a

沙哑音 鼻音

个人特点 声学特点 语音习惯 沙哑音 1

1/5

1/5

1/9

1/7

鼻音 1/5 1 1 1/5 1/5 个人特点 1/5

1

1

1/5

1/5

声学特点 1/9 1/5 1/5 1 3

语音习惯

1/7 1/5 1/5 3 1

从而得到正反矩阵：

A=???

????

?????

?

??

?13/1557

315595/15/11155/15/11157/19/15/15/11

利用MATLAB 编程求得矩阵A 的最大特征值;

对正互反矩阵进行一致性检验，采用T.L.Saaty 一致性指标：CI=1

--n n

λ，

一致性对比率CR=

RI CI =12

.1087.0=0.0777<0.1,即通过一致性检验。 标准化特征向量a (2)=(0.0308 0.0927 0.0927 0.4839 0.2999)

即当前沙哑音，鼻音，个人特点，声学特征，语音习惯的权重为：

(0.0308 0.0927 0.0927 0.4839 0.2999)

下面构造方案层对准则层的每个准则的正互反矩阵：

B 1=?????

??

??

?

??

???

??

???11555

71155575/15/1111

35/15/111135/15/111137/17/13/13/13/11 B 2=?????

?

?

?????13/13/15/13113/13113/15331 B 4=??????????173/17/115/13/151 B 5=??

??

?

?????13/113/11313/11 所有CR k 均小于0.1，均通过一致性检验。

我们把各方案对目标的权向量，称为组合权向量。记做

W (2)=??

?

??

??

?

?

???

?

?????

??

?

??

???005.00781.0025.05017.05.01998.005.00953.001998.03598.025.04030.005222.03598.000000814.000000814.000000814.000000361

.0 所以方案C k 在目标中的组合权中应该为a (2)与a (3)对应向量的两两乘积之

和，因此可得C 1 ,C 2 ,C 3 ,C 4 ,C 5 ,C 6 ,C 7 ,C 8在目标中的组合权重。从而得到方案层对准则层组合权向量：

a (3)= (0.0011,0.0025,0.0025,0.3295,0.2257,0.3826,0.0536)

结论：

由当前权重比例来看，建立代表说话人身份的语音特征模型中的8个影响因素声学特征占有最大的权重。而语者的声学特征能用Mel 倒频谱系数(MFCC)很好的描述。

所以接下来通过matlab 编程从源语音数据中提取说话人特征（Mel 倒频谱系数）从而建立代表语者身份的语音特征模型。 4.1.2 语音特征模型的建立与求解

本文采用识别性能较好且能充分模拟人耳听觉感知特性的Mel 频率倒谱系数（MFCC ）,提取的各过程如下：

图1

MFCC 系数提取流程图

1. 由于语音信号的时变特性使得特征参数提取必须在一小段语音上进行，

因此在特征参数提取前先对采集到的语音信号进行预加重、加窗和分帧等预处理。此部分的处理用matlab 完成，程序代码在附录中。

2. 将预处理后的语音信号x(n)进行快速傅里叶转变(FFT)得到频谱X(k):

()N-1

2/0

X k x(),01j k N n n e k N π-==≤≤-∑ (1)

其中，N 为傅里叶变换的样点数。本论文中N 取256

3. 求频谱X(k)的平方，即能量谱，然后通过Mel 滤波器H q (k)对频谱进行

平滑，并消除谐波，凸显原先语音的共振峰。H q (k)是一组三角带通滤波器，设计如下：

在本论文中区滤波器个数为Q=24，它的图像如下：

图2

Mel 频率滤波器组

()q 0,(1)(1),(1)()()(1)

H k (1),()(1)(1)()0,(1)

k f q k f q f q k f q f q f q f q k f q k f q f q f q k f q <-??--?-≤<--?=?

+-?≤≤+?+-?

>+? （2）

式中：f(q),q=1,…,Q 是中心频率，Q 是滤波器个数。

4. 对滤波器组的输出取对数：压缩语音频谱的动态范围；将频域中噪声的

乘性成分转换为加性成分，对数Mel 频谱S(q)为：

()()()210S q ln N q k X k H k -=????

=??????

∑

（3）

5. 离散余弦变换（DCT ）

将（3）获得的S(q)变换到时域，结果就是Mel 频率倒谱系数(MFCC),第n 个系数C(n)的计算如下式：

()()()Q 1

q 0n q 0.52C n S q cos 0n

∑， （4） 其中，L 是MFCC 阶数，一般为12~16，Q 为23~26.

MFCC 反映了语音信号的静态特性，随时间的变化反映了语音特性的动态特性，用差分公式表示如下：

()2

m C n ()M M

m i M i M

C n i i i =-=-??

?=+??

??

∑∑ （5）

这里M 的值一般取2或3。

至此，MFCC 特征提取结束。在实际过程中，可以根据需要与实验测试来选取维数不同的特征向量，如：12，13，24，26，36，39维等。本文选用24维的MFCC 特征向量来作为声纹识别系统的特征提取结果。 该模型的matlab 代码见附录一。 4.2 问题二模型的建立和求解：

4.2.1基于GMM 的声纹识别建立的模型

GMM 本质上是利用多维概率密度函数对语言信号进行建模。一个具有M 阶混合分量的D 维GMM 表示如下

()()1

P X M

i i i w f X λ==∑

（6）

式(6)表明了样本X 属于类别λ的概率大小，M 实际上是表示样本X 的类别数。 式(6)中，i w 为混合权重，且1

1M

i

i w ==∑。()i

f X 是D 维的联合高斯分布，表

示为

()()()()1T 1exp 22i i i D i

X u X u f X π-??

---??

??=

∑∑

（7）

式中，i u 是均值，

i

∑

是协方差矩阵。GMM 的所有参数用λ 表述为

{},, 1,,i i i u w i M λ==∑ （8）

4.2.2

GMM 模型的求解

4.2.2.1 初始值

由k 均值(k-mean)聚类算法对样本进行聚类，利用各类的均值作为0i u ，

并计算

i ∑

，0i w 取各类样本占样本总数的比例。

4.2.2.2 最大化步骤

对于一组长度为T 的语言序列{}12X ,,,T X X X = ，它的似然度函数表示

为

()()1

P P T

t t X X λλ==∏

（9）

为说话人建立GMM 后，我们用EM 算法估计GMM 的所有参数，通过迭代

找到一个新的模型*

λ ，使()()*P X P X λλ≥。最终得到GMM 各参数的

重估公式： 更新混合权重

()1

1,T

i t t w P i X T λ==∑

（10）

● 更新均值

()()

11

,,T

t t

t i T

t t P i

X X u P i

X λλ===

∑∑

（11）

● 更新方差矩阵

()()()

2

1

1

,,T

t t i

t i

T

t t P i X X u P i

X λλ==-=

∑∑

∑ （12）

4.2.2.3 估计步骤

对基于GMM 的说话人识别，采用帧似然概率得分法进行判决。具有N

各说话人的识别系统，相应的GMM 分别为12,,,N λλλ 。在识别时，给定

待识别语音序列{}12X

,,,T X X X = ，则这段语音属于第n 个说话人的最

大后验概率为：

()(|)()

P ()

n n n P X P X P X λλλ=

（13）

式中，P(X)为所有说话人条件下X 的概率密度，对于确定的X ，P(X)是一个确定

的常数，()P

n λ 为第n 个人说话的先验概率，且

()1

P ,1n n N N

λ=≤≤ （14）

4.2.2.4 收敛条件

不断地迭代2.2.2和2.2.3步骤，重复更新(10) (11) (12)三个式子的值，直到

()*()P X P X λλε-<

（15）

式(15)中的ε通常取10-5，()P X λ通过公式(9)计算，*()P X λ表示更新参

数的值。

4.2.2.5 识别结果的得出

()P X|n λ 为X 是第n 个人说话的条件概率，因此求取后概率的最大值可

以通过求取()P

X|n λ获得，识别结果为：

()*max |n n P X λ=

（16）

4.3

4.4 问题四的解答：

我们知道声纹识别需要解决的关键问题还有很多，诸如：短话音问题，能否用很短的语音进行模型训练，而且用很短的时间进行识别，这主要是声音不易获取的应用所需求的；声音模仿(或放录音)问题，要有效地区分开模仿声音(录音)和真正的声音；多说话人情况下目标说话人的有效检出；消除或减弱声音变化(不同语言、内容、方式、身体状况、时间、年龄等)带来的影响；消除信道差异和背景噪音带来的影响;此时需要用到其他一些技术来辅助完成，如去噪、自适应等技术。

设长度为N 的M 维特征矢量序列为：X={X 1 ，X 2 ，?，X N )，其中第n(0

2

2

1

1S (X

)

M

nm

n m X M

==

-∑ （17）

式中： 为特征矢量的维数； 为第，2个矢量的第m 个参数； 为第n 个矢量的平均值，第m 维矢量的权值为:

2

1m

m S =π （18）

相应的基于方差的加权欧氏距离公式为:

∑--

=

M

m

km

nm

m

n

C

X

K

X

D

1

2

) (

)

,

(π（19）

式中：C km为第K个类的聚类中心。

方差是用来衡量样本数据的离散程度的，进行聚类分析时易受样本孤立点的响。为了消除此影响，文中利用距离和对样本的孤立点进行了处理计算矢量集中矢量两两之间的几何距离，累加距离和，去掉L个与其他矢量之问距离和最大的矢量，然后从剩余的矢量中选出距离最大的2个矢量作为不同类的聚类中心，着从其余的矢量中寻找出到所有聚类中心的距离和最大的矢量作为另一个类的中心，直到选出K个聚类中心。改进后的K—means聚类算法描述如下：

(1)从上述方法选出的K个聚类中心出发，利用公式(19)计算矢量集中的矢量与各类中心的距离，把剩余矢量划分到离它距离最近的类中，形成初始聚类；

(2)按照步骤(1)的聚类，更新各个类的中心；

(3)以新聚类中心为参照点迭代步骤(1)和(2)，到聚类中心不再变化或变化微小时停止。

五、结论分析与推广

声纹识别可以应用的范围很宽，可以说声纹识别几乎可以应用到人们日常生活的各个角落。比如下面举几个例子。

(1)信息领域。比如在自动总机系统中，把身份证之声纹辨认和关键词检出器结合起来，可以在姓名自动拨号的同时向受话方提供主叫方的身份信息。前者用于身份认证，后者用于内容认证。同样，声纹识别技术可以在呼叫中心(Call Center)应用中为注册的常客户提供友好的个性化服务。

(2)银行、证券。鉴于密码的安全性不高，可以用声纹识别技术对电话银行、远程炒股等业务中的用户身份进行确认，为了提供安全性，还可以采取一些其他措施，如密码和声纹双保险，如随机提示文本用文本相关的声纹识别技术进行身份确认(随机提示文本保证无法用事先录好的音去假冒)，甚至可以把交易时的声音录下来以备查询。

(3)公安司法。对于各种电话勒索、绑架、电话人身攻击等案件，声纹辨认技术可以在一段录音中查找出嫌疑人或缩小侦察范围；声

纹确认技术还可以在法庭上提供身份确认（同一性鉴定）的旁证。在监狱亲情电话应用中，通过采集犯人家属的声纹信息，可有效鉴别家属身份的合法性。在司法社区矫正应用中，通过识别定位手机位置和呼叫对象说话声音的个人特征，系统就可以快速的自动判断被监控人是否在规定的时间出现在规定的场所，有效地解决人机分离问题。

(4)军队和国防。声纹辨认技术可以察觉电话交谈过程中是否有关键说话人出现，继而对交谈的内容进行跟踪(战场环境监听)；在通过电话发出军事指令时，可以对发出命令的人的身份进行确认(敌我指战员鉴别)。目前该技术在国外军事方面已经有所应用，据报道，迫降在我国海南机场的美军EP-3侦察机中就载有类似的声纹识别侦听模块。

(5)保安和证件防伪。如机密场所的门禁系统。又如声纹识别确认可用于信用卡、银行自动取款机、门、车的钥匙卡、授权使用的电脑、声纹锁以及特殊通道口的身份卡，把声纹存在卡上，在需要时，持卡者只要将卡插入专用机的插口上，通过一个传声器读出事先已储存的暗码，同时仪器接收持卡者发出的声音，然后进行分析比较，从而完成身份确认。同样可以把含有某人声纹特征的芯片嵌入到证件之中，通过上面所述的过程完成证件防伪。

介于应用方面的考虑，该模型的建立对声纹识别技术在信息领域，银行、证券，公安司法，军队和国防，保安和证件防伪，社保基金等方面的正常发展和实施具有里程碑的的作用，也为以后的有关身份证实问题的行业的进步具有重要意义。除此以外声纹识别所提供的安全性可与其他生物识别技术（如：指纹、掌形和虹膜）相媲美，而且语音采集装置造价低廉，只需电话/手机或麦克风即可，无需特殊的设备；它与说话语言无关，与方言腔调无关，不涉及隐私问题，适应人群范围很广；声音信号便于远程传输和获取，在基于电信和网络的身份识别应用中，声纹识别更为擅长，得天独厚。

六、模型评价

优点：MFCC模型应用广泛，是对人耳的仿真，准确度高。GMM模型能很好的仿真语者广阔的语音空间中的某一类下的一项声学类别，对任一种具有多类别的样本而言，它都具有极佳的能力去近似这些样本（在声纹识别中即为特征向量）空间概率分布。

不足：传统的K-means算法在聚类过程中采取几何距离就近原则，这只有当类的自然分布为球状或接近于球状时，也就是说当类中各矢量的分布接近于相等时才能取得较好的识别效果。而语音信号是近似椭圆状的正态分布，各矢量的分布不等，若用传统K-means算法作为GMM 的初始化方法，将会忽略特征矢量各维分量对识别的影响，无法得到令人满，意的识别结果。

改进：(1)从上述方法选出的个聚类中心出发，计算矢量集中的矢量与各类中心的距离，把剩余矢量划分到离它距离最近的类中，形成初始聚类；

(2)按照步骤(1)的聚类，更新各个类的中心；

(3)以新聚类中心为参照点迭代步骤(1)和(2)，直到聚类中心不再变化

或变化微小时停止。

模型的推广：特征参数提取方面，本文在语音信号特征提取方面仅

选用了语音领域较为认同的美尔倒频谱系数（MFCC），而目前小波变换

以其时域、频域兼而有之的高分辨率优势正在不断被引进到语音信号的

特征提取过程中。相信进一步的研究，必定会得到不少的收获。本文并

没有特别深入地针对某些比较细节的问题来做处理，而在实际的声纹识

别系统中，有许多世纪的因素会影响系统的表现，例如语者模型老化、

恶劣的声学环境等，这些问题都有利于进一步的深入研究。

七、参考文献

[1]冷建华，李萍，王良红数字信号处理北京：国防工业出版社，2002.

[2]滕素珍数理统计大连：大连理工大学出版社，2008.

[3]吴金池语者辨识系统之研究国立中央大学电机工程研究所硕士论文，2002。

[4]卢文祥，杜润生机械工程测试·信号·信号分析武汉：华中科技大学出版社，1999.

[5]易克初，田斌，付强语言信号处理北京：国防工业出版社，2000.

[6]石燕声纹识别技术研究南京航空航天大学硕士学位论文，2004.

[7]朱民雄计算机语言技术[M] 北京：北京航空航天大学出版社，1992.

[8]Juang B.语音处理的过去、现在和未来[J]，IEEE信号处理，1998.

[9]吴玺洪声纹识别听声辨认[N] 计算机世界，2001.

[10]张彩娟，霍春宝改进K-means算法在声纹识别中的应用，2011年10月

附录1:

%---------------demo1.m-------------------------%

x=wavread('001.wav');%读入源语音文件

mel_data=mfcc(x); %输出mel频率倒频系数

%--------------------函数mfcc.m的源代码--------------------%

function ccc=mfcc(x)

%% x为用wavread函数读到的初始信号值

%归一化mel滤波器组系数

bank=melbankm(24,256,8000,0,0.5,'m');

bank=full(bank);

bank=bank/max(bank(:));

%DCT系数，12*24

for k=1:12

n=0:23;

dctcoef(k,:)=cos((2*n+1)*k*pi/(2*24));

end

%归一化倒谱提升窗口

w=1+6*sin(pi*[1:12]./12);

w=w/max(w);

%预加重滤波器

xx=double(x);

xx=filter([1 -1],1,xx);

%语音信号分帧

%xx=enframe(xx,256,80);

xppl=length(xx);

j=1;

for i=65:80:xppl-256,

xx1(j,:)=xx(i:i+256-1)';

j=j+1;

end

xx=xx1;

%计算每帧的MFCC参数

for i=1:size(xx,1)

y=xx(i,:);

s=y'.*hamming(256);

t=abs(fft(s));

t=t.^2;

t=t+2*realmin;

c1=dctcoef*log(bank*t(1:129));

c2=c1.*w';

m(i,:)=c2';

end

%差分参数

dtm=zeros(size(m));

for i=3:size(m,1)-2

dtm(i,:)=-2*m(i-2,:)-m(i-1,:)+m(i+1,:)+2*m(i+2,:);

end

dtm=dtm/3;

%合并mfcc参数和一阶差分mfcc参数

ccc=[m dtm];

%去除首尾两帧，因为这两帧的一阶差分参数为0

ccc=ccc(3:size(m,1)-2,:);

return

%%----------------enframe函数的matlab代码----------------------%%% function [f,t,w]=enframe(x,win,inc,m)

%ENFRAME split signal up into (overlapping) frames: one per row. [F,T]=(X,WIN,INC)

%

% Usage: (1) f=enframe(x,n) % split into frames of length n %

% (2) f=enframe(x,hamming(n,'periodic'),n/4) % use a 75% overlapped Hamming window of length n

%

% Inputs: x input signal

% win window or window length in samples

% inc frame increment in samples

% m mode input:

% 'z' zero pad to fill up final frame

% 'r' reflect last few samples for final frame % 'A' calculate window times as the centre of mass % 'E' calculate window times as the centre of energy %

% Outputs: f enframed data - one frame per row

% t fractional time in samples at the centre of each frame % w window function used

%

% By default, the number of frames will be rounded down to the nearest % integer and the last few samples of x() will be ignored unless its length % is lw more than a multiple of inc. If the 'z' or 'r' options are given, % the number of frame will instead be rounded up and no samples will be ignored.

%

% Example of frame-based processing:

% INC=20 % set frame increment in samples

% NW=INC*2 % oversample by a factor of 2 (4 is also often used)

% S=cos((0:NW*7)*6*pi/NW); % example input signal

% W=sqrt(hamming(NW),'periodic')); % sqrt hamming window of period NW

% F=enframe(S,W,INC); % split into frames % ... process frames ...

% X=overlapadd(F,W,INC); % reconstitute the time waveform (omit "X=" to plot waveform)

% Bugs/Suggestions:

% (1) Possible additional mode options:

% 'u' modify window for first and last few frames to ensure WOLA % 'a' normalize window to give a mean of unity after overlaps % 'e' normalize window to give an energy of unity after overlaps

% 'wm' use Hamming window

% 'wn' use Hanning window

% 'x' include all frames that include any of the x samples

% Copyright (C) Mike Brookes 1997-2012

% Version: $Id: enframe.m 1713 2012-03-30 21:27:46Z dmb $

%

% VOICEBOX is a MATLAB toolbox for speech processing.

% Home page: https://www.360docs.net/doc/2f12565677.html,/hp/staff/dmb/voicebox/voicebox.html

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%

% This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by

% the Free Software Foundation; either version 2 of the License, or % (at your option) any later version.

%

% This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details.

%

% You can obtain a copy of the GNU General Public License from % https://www.360docs.net/doc/2f12565677.html,/copyleft/gpl.html or by writing to

% Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%

nx=length(x(:));

if nargin<2 || isempty(win)

win=nx;

end

if nargin<4 || isempty(m)

m='';

end

nwin=length(win);

if nwin == 1

lw = win;

w = ones(1,lw);

else

lw = nwin;

w = win(:)';

end

if (nargin < 3) || isempty(inc)

inc = lw;

end

nli=nx-lw+inc;

nf = fix((nli)/inc);

na=nli-inc*nf;

f=zeros(nf,lw);

indf= inc*(0:(nf-1)).';

inds = (1:lw);

f(:) = x(indf(:,ones(1,lw))+inds(ones(nf,1),:));

if nargin>3 && (any(m=='z') || any(m=='r')) && na>0

if any(m=='r')

ix=1+mod(nx-na:nx-na+lw-1,2*nx);

f(nf+1,:)=x(ix+(ix>nx).*(2*nx+1-2*ix));

else

f(nf+1,1:na)=x(1+nx-na:nx);

end

nf=size(f,1);

end

if (nwin > 1) % if we have a non-unity window

f = f .* w(ones(nf,1),:);

end

if nargout>1

if any(m=='E')

t0=sum((1:lw).*w.^2)/sum(w.^2);

elseif any(m=='E')

t0=sum((1:lw).*w)/sum(w);

else

t0=(1+lw)/2;

end

t=t0+inc*(0:(nf-1)).';

end

%%%-----------melbankm函数的代码-----------%%

function [x,mn,mx]=melbankm(p,n,fs,fl,fh,w)

%MELBANKM determine matrix for a mel-spaced filterbank [X,MN,MX]=(P,N,FS,FL,FH,W)

%

% Inputs: p number of filters in filterbank

% n length of fft

% fs sample rate in Hz

% fl low end of the lowest filter as a fraction of fs (default = 0)

% fh high end of highest filter as a fraction of fs (default =

0.5)

% w any sensible combination of the following:

% 't' triangular shaped filters in mel domain (default)

% 'n' hanning shaped filters in mel domain

% 'm' hamming shaped filters in mel domain

%

% 'z' highest and lowest filters taper down to zero (default)

% 'y' lowest filter remains at 1 down to 0 frequency and

% highest filter remains at 1 up to nyquist freqency

%

% If 'ty' or 'ny' is specified, the total power in the fft

is preserved.

%

% Outputs: x a sparse matrix containing the filterbank amplitudes

% If x is the only output argument then

size(x)=[p,1+floor(n/2)]

% otherwise size(x)=[p,mx-mn+1]

% mn the lowest fft bin with a non-zero coefficient

% mx the highest fft bin with a non-zero coefficient

%

% Usage: f=fft(s); f=fft(s);

% x=melbankm(p,n,fs); [x,na,nb]=melbankm(p,n,fs);

% n2=1+floor(n/2); z=log(x*(f(na:nb)).*conj(f(na:nb)));

% z=log(x*abs(f(1:n2)).^2);

% c=dct(z); c(1)=[];

%

% To plot filterbanks e.g. plot(melbankm(20,256,8000)')

%

% Copyright (C) Mike Brookes 1997

% Version: $Id: melbankm.m,v 1.3 2005/02/21 15:22:13 dmb Exp $

%

% VOICEBOX is a MATLAB toolbox for speech processing.

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%

% This program is free software; you can redistribute it and/or modify

% it under the terms of the GNU General Public License as published

by

% the Free Software Foundation; either version 2 of the License, or % (at your option) any later version.

%

% This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details.

%

% You can obtain a copy of the GNU General Public License from % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%

if nargin < 6

w='tz';

if nargin < 5

fh=0.5;

if nargin < 4

fl=0;

end

end

end

f0=700/fs;

fn2=floor(n/2);

lr=log((f0+fh)/(f0+fl))/(p+1);

% convert to fft bin numbers with 0 for DC term

bl=n*((f0+fl)*exp([0 1 p p+1]*lr)-f0);

b2=ceil(bl(2));

b3=floor(bl(3));

if any(w=='y')

pf=log((f0+(b2:b3)/n)/(f0+fl))/lr;

fp=floor(pf);

r=[ones(1,b2) fp fp+1 p*ones(1,fn2-b3)];

c=[1:b3+1 b2+1:fn2+1];

v=2*[0.5 ones(1,b2-1) 1-pf+fp pf-fp ones(1,fn2-b3-1) 0.5];

mn=1;

mx=fn2+1;

else

b1=floor(bl(1))+1;

b4=min(fn2,ceil(bl(4)))-1;

pf=log((f0+(b1:b4)/n)/(f0+fl))/lr;

fp=floor(pf);

pm=pf-fp;