当前位置:首页 >> 工学 >> 数字语音信号处理实验(学生)

数字语音信号处理实验(学生)


数字语音信号处理 实验指导书

北方学院信息科学与工程学院 电子教研室

2014 年 1 月

前言

语音信号处理是研究用数字信号处理技术和语音学知识对语音信号进行处理的新兴的学科,是目前发展最 为迅速的信息科学研究领域的核心技术之一。通过语音传递信息是人类最重要、最有效、最常用和最方便的交 换信息形式。同时,语言也是人与机器之间进行通信的重要工具,它是一种理想的人机通信方式,因而可为信 息处理系统建立良好的人机交互环境,进一步推动计算机和其他智能机器的应用,提高社会的信息化程度。 语音信号处理是一门新兴的学科,同时又是综合性的多学科领域和涉及面很广的交叉学科。虽然从事这一 领域研究的人员主要来自信号与信息处理及计算机应用等学科,但是它与语音学、语言学、声学、认知科学、 生理学、心理学等许多学科也有非常密切的联系。 20 世纪 60 年代中期形成的一系列数字信号处理的理论和算法,如数字滤波器、快速傅立叶变换(FFT)等 是语音信号数字处理的理论和技术基础。随着信息科学技术的飞速发展,语音信号处理取得了重大的进展:进 入 70 年代之后,提出了用于语音信号的信息压缩和特征提取的线性预测技术(LPC) ,并已成为语音信号处理 最强有力的工具,广泛应用于语音信号的分析、合成及各个应用领域,以及用于输入语音与参考样本之间时间 匹配的动态规划方法;80 年代初一种新的基于聚类分析的高效数据压缩技术—矢量量化(VQ)应用于语音信 号处理中;而用隐马尔可夫模型(HMM)描述语音信号过程的产生是 80 年代语音信号处理技术的重大发展, 目前 HMM 已构成了现代语音识别研究的重要基石。近年来人工神经网络(ANN)的研究取得了迅速发展,语音 信号处理的各项课题是促进其发展的重要动力之一,同时,它的许多成果也体现在有关语音信号处理的各项技 术之中。 为了深入理解语音信号数字处理的基础理论、算法原理、研究方法和难点,根据数字语音信号处理教学大 纲,结合课程建设的需求,我们编写了本实验参考书。 本参考书针对教学大纲规定的八个研究设计型实验,每个实验给出了参考程序,目的是起一个抛砖引玉的 作用,学生在学习过程中,可以针对某一个实验进行延伸的创新学习,比如说,语音端点的检测、语音共振峰 提取、基于 HMM 或 DTW 的有限词汇或大词汇的特定人、非特定人的语音识别、识别率的提高(如何提高有 噪环境下的识别率) 、以及编码问题等,同时在学习中还可深入思考如何将有关的方法在嵌入式系统或 DSP 下 的实现问题等。

第一章 数字语音信号处理教学大纲

一、课程说明
学 分 数:2.5 总 学 时:54 学时分配:讲课 36 学时,实验 18 学时 适用专业:电子信息科学与技术

二、课程教学的目的与任务
本课程的学习目的是掌握语音信号处理的基本理论、基本分析方法,了解在语音信号处理领域中相关研究 热点,激发学习者对语音处理相关研究方向中的有关兴趣,为以后的开展语音处理相关领域的研究、开发打下 一个良好的基础。本课程是电子信息科学与技术专业的方向模块课。 本门课程的教学分理论和实验教学两部分,理论教学注重培养学生基本问题的分析方法,从而掌握基本的 语音信号处理的理论与概念,理论教学还包括多种形式的自主学习,如网上学习、课外阅读、大型作业、主题 调查、读书报告、分组讨论等。实验教学注重培养学生的动手能力、分析和解决问题的能力。

三、课程教学的基本内容及学时分配
1. 语音信号处理概述 语音信号处理的发展概况,语音信号处理的应用。 2.语音信号的特性及模型 (理论教学:2 学时) (理论教学:2 学时)

语音信号的产生,语音信号的特性,语音信号产生的数字模型,语音感知。 3. 语音信号的时域分析 (理论教学:2 学时,自主学习:2 学时)

语音信号的数字化和预处理,短时能量分析,短时过零分析,短时相关分析。 4. 语音信号的频域分析 (理论教学:2 学时,自主学习:2 学时)

短时傅里叶变换,短时傅里叶变换的取样率,语音信号的短时综合,语谱图。 5. 语音信号的同态滤波及倒谱分析 (理论教学:2 学时,自主学习:1 学时)

同态信号处理的基本原理,复倒谱和倒谱,语音信号两个卷积分量复倒谱的性质,避免相位卷绕的算法, 语音信号复倒谱分析实例。 6. 语音信号的 LP 分析 (理论教学:2 学时,自主学习:2 学时)

线性预测分析的基本原理, 线性预测方程组的建立, 线性预测分析的解法, 线性预测分析应用, 线谱对(LSP) 分析,极零模型。 7. 语音信号的矢量量化 (理论教学:2 学时,自主学习:1 学时)

矢量量化的基本原理,失真测度,最佳矢量量化器和码本的设计,降低复杂度的矢量量化系统,语音参数

的矢量量化。 8. 语音编码-波形编码法 (理论教学:2 学时,自主学习:0.5 学时)

语音信号的压缩编码原理,脉冲编码调制(PCM)及其自适应,预测编码及其自适应 APC,自适应差分脉冲 编码调制(ADPCM)及自适应增量调制(ADM),子带编码(SBC),自适应变换编码(ATC)。 9. 语音编码-参数编码法 (理论教学:2 学时,自主学习:0.5 学时)

声码器的基本结构,相位声码器和通道声码器,同态声码器, 线性预测声码器,混合编码,各种语音编码 方法的比较,语音编码的性能指标和质量评价。 10. 隐马尔可夫模型(HMM) (理论教学:2 学时,自主学习:1 学时)

隐马尔可夫模型的引入,隐马尔可夫模型的定义,隐马尔可夫模型三项问题的求解,HMM 的一些实际问 题。 11. 语音识别技术 (理论教学:2 学时,自主学习:1 学时)

语音识别概述,动态时间规整(DTW)识别技术,隐马尔可夫模型(HMM)识别技术,语音识别的应用 技术。 12. 语音合成、语音增强技术 (理论教学:2 学时,自主学习:1 学时)

语音合成原理,共振峰合成,线性预测合成,专用语音合成硬件及语音合成器芯片,语音增强。

四、教学方法
本课程总学时 54(总学分:2.5) ;其中课堂讲授:24 学时;自主学习:12 学时;实验:18 学时。 理论课采用课堂教学方式,使用多媒体辅助教学手段,进行基本内容的讲授。适当安排一定的习题课时间, 并布置适当的设计题以培养学生的设计、分析问题的能力。 自主学习内容由学生自主学习参考教材的内容,并采用多种渠道,如查阅最新语音信号处理方面的科技文 献、资料,作出学习报告。目的是培养学生的自学能力和科技文献的检索和查阅能力,同时可以有助于学生了 解和掌握语音信号处理领域的最新技术进展和应用情况,将理论知识和实际应用结合起来,促进学生学习的积 极性和主动性。本课程讲授自主学习的内容依每部分的教学进度交替安排。 实验为研究型(设计型)实验,共安排 8 个实验,为了真正达到研究设计型实验的目的,将自主学习和研 究设计型实验结合起来,统一安排。

五、考核及成绩评定方式
本课程的考核内容由下面四部分组成: 1、期末考试 M1(100 分) 考核内容:教学计划全部内容;考核形式:闭卷考试。占总评成绩的 70% 2、实验考核(含自主学习)M2(100 分) 设计型实验共占 10%,评分标准是按试验分析方法、所设计的实验程序、实验结果等,由任课教师评定成

绩 3、论文及主题报告 M3(100 分) 按一般科学论文的写作规范的要求,写作专题论文(含自主学习) ,每一学生选择至少一个写作规范的专题 论文进行课堂交流报告,根据论文写作水平、报告的内容、思路、对问题的理解、以及报告方式等评定成绩。 4、平时考核 M4(100 分) 由任课主讲教师按课堂表现、平时实验、自主学习情况及作业评定成绩 期末总评成绩 M=M1× 60%+M2× 10%+M3× 20%+M4× 10%

六、教材及参考书目
推荐教材:张雄伟等编著, 《现代语音处理技术及应用》 ,机械工业出版社,2003 年。 参考教材: 1、 L.R. Rabiner, B.H. Juang. Fundamentals of Speech Recognition. Prentice Hall, Englewood Cliffs,1993. 清华 大学出版社(影印) ,2002 年. 2、 胡航. 语音信号处理(修订版) ,哈尔滨工业大学出版社,2002 年. 3、 易克初,田 斌等. 语音信号处理,国防工业出版社,2000 年. 4、 赵 力. 语音信号处理,机械工业出版社,2003 年. 5、 吴家安等. 语音编码技术及应用,机械工业出版社,2006 年. 6、 韩继庆,张 磊,郑铁然. 语音信号处理,清华大学出版社,2004 年. 7、 D.G.Childers. Matlab 之语音处理与合成工具箱(影印版) ,清华大学出版社,2004 年. 8、 Thomas F. Quatieri 著,赵胜辉等译, 《离散时间语音信号处理—原理与应用》 ,电子工业出版社,2004.

七、实践环节
实验学时数:18 实验项目数:8 1、目的与基本要求 实验为研究型(设计型)实验,共安排 8 个,为了真正达到研究设计型实验的目的,采用开放实验的办法, 将自主学习和研究设计型实验结合起来,统一安排。 通过开放实验,目的使学生进一步理解数字语音信息处理的基本方法,提高学生自主分析、发现及解决问 题的能力,锻炼学生论文写作能力,为实际的应用打下扎实的基础。 2、研究设计型实验的内容 1)研究设计型实验 1: 语音信号的采集及预处理 2)研究设计型实验 2: 语音及音乐信号的采样、滤波

3)研究设计型实验 3 基于 MATLAB 的语音信号时域特征分析 要求: 按所学相关语音处理得的知识,通过网上学习、资料查阅,自己设计程序,给出某一语音信号的短时过零 率、短时能量、短时自相关特征的分析结果,并借助时域分析方法检测所分析语音信号的基音周期,写出报告 (按一般科学论文的写作规范) 。 4)研究设计型实验 4: 基于 MATLAB 分析语音信号频域特征 要求: 按所学相关语音处理的得知识,通过网上学习、资料查阅,自己设计程序,给出某一语音信号的短时谱、 倒谱、语谱图的分析结果,并借助频域分析方法检测所分析语音信号的基音周期或共振峰,写出报告(按一般 科学论文的写作规范) 。 5)研究设计型实验 5: 基于 MATLAB 进行语音信号的 LPC 分析 要求: 按所学相关语音处理的知识,通过网上学习、资料查阅,自己设计程序,给出某一语音信号的 LPC 分析结 果,包括 LPC 谱、LPCC 谱的分析结果,并借助 LPC 分析方法检测所分析语音信号的基音周期和共振峰,写出 报告(按一般科学论文的写作规范) 。 6)研究设计型实验 6: 基于 VQ 的特定人孤立词语音识别研究 要求: 按所学相关语音处理的知识,通过网上学习、资料查阅,借助 MATLAB 工具,自己设计基于 VQ 的码本训 练程序和识别程序(尽量选用所学 HMM 或 DTW 方法设计识别程序) ,能识别特定人的语音,分析所设计系统 的特性,写出报告(按一般科学论文的写作规范) 。 7)研究设计型实验 7: 音乐信号处理 8)研究设计型实验 8: 双音多频(DTMF)信号的合成和识别

第二章

实验

实验一 语音信号的采集及预处理
一、实验目的 在理论学习的基础上, 进一步地理解和掌握语音信号预处理及短时加窗的意义及基于 matlab 的实现方法。 二、实验原理 1. 语音信号的录音、读入、放音等:练习 matlab 中几个音频处理函数,利用函数 wavread 对语音信号进行采 样,记住采样频率和采样点数,给出以下语音的波形图(2.wav) 。利用 wavplay 或 soundview 放音。也可以 利用 wavrecord 自己录制一段语音,并进行以上操作(需要话筒)。 2. 语音信号的分帧:对语音信号进行分帧,可以利用 voicebox 工具箱中的函数 enframe。voicebox 工具箱是 基于 GNU 协议的自由软件,其中包含了很多语音信号相关的函数。 3.语音信号的加窗: 本步要求利用 window 函数设计窗口长度为 256(N=256)的矩形窗(rectwin)、 汉明窗(hamming) 及汉宁窗(hann)),利用 wvtool 函数观察其时域波形图及频谱特性,比较得出结论。观察整个信号加矩形窗 及汉明窗后的波形,利用 subplot 与 reshape 函数将分帧后波形、加矩形窗波形及加汉明窗波形画在一张图 上比较。取出其中一帧,利用 subplot 与 reshape 函数将一帧语音的波形、加矩形窗波形及加汉明窗波形画 在一张图上比较将得出结论。

z 。 4. 预加重:即语音信号通过一个一阶高通滤波器 1 ? 0.9375
三、实验步骤、实验程序、图形及结论 1.语音信号的录音、读入、放音等 程序: [x,fs,nbit]=wavread('D:\2.wav'); %fs=10000,nbit=16 y=soundview('D:\2.wav') 2.语音信号的分帧 程序: [x,fs,nbit]=wavread('D:\2.wav'); len=256; inc=128; y=enframe(x,len,inc); figure; subplot(2,1,1),plot(x) subplot(2,1,2),plot(y)

?1

3.语音信号加窗: 程序: N=256; w = window('rectangle',N); w1 = window('hamming',N); w2 = window('hanning',N); wvtool(w,w1,w2)

4.预加重 程序: [x,fs,nbit]=wavread('D:\2.wav'); len=256; inc=128; y=enframe(x,len,inc); z=filter([1-0.9375],1,y) figure(2) subplot(2,1,1),plot(y) subplot(2,1,2),plot(z)

四、思考题 1.语音信号包括哪些预处理,作用分别是什么? 2.不同窗口的优缺点,窗口长度如何选取?

实验二
1、实验目的

语音及音乐信号的采样、滤波

(1) 理解采样率和量化级数对语音信号的影响; (2) 设计滤波器解决实际问题。 2、 实验要求 利用电脑的声卡录一段语音信号及音乐信号, (1) 观察使用不同采样率及量化级数所得到的信号的听觉效果,从而确定对不同信号的最佳的 采样率; (2) 分析音乐信号的采样率为什么要比语音的采样率高才能得到较好的听觉效果; (3) 注意观察信号中的噪声(特别是 50hz 交流电信号对录音的干扰,设计一个滤波器去除该噪 声。 3、 实验提示 (1) 推荐录音及播放软件:CoolEdit;

(2) 分析语音及音乐信号的频谱,根据信号的频率特性理解采样定律对信号数字化的工程指导

意义; (3) 可用带阻滤波器对 50Hz 交流电噪声进行去噪处理; (4) 也可研究设计自适应滤波器对 50Hz 噪声及其它随机环境噪声进行滤波处理。

实验三 基于 MATLAB 的语音信号时域特征分析
一、实验目的
语音信号是一种非平稳的时变信号,它携带着各种信息。在语音编码、语音合成、语音识别和语音增强等 语音处理中无一例外需要提取语音中包含的各种信息。 语音信号分析的目的就在与方便有效的提取并表示语音 信号所携带的信息。语音信号分析可以分为时域和变换域等处理方法,其中时域分析是最简单的方法,直接对 语音信号的时域波形进行分析, 提取的特征参数主要有语音的短时能量, 短时平均过零率, 短时自相关函数等。 本实验要求掌握时域特征分析原理,并利用已学知识,编写程序求解语音信号的短时过零率、短时能量、 短时自相关特征,分析实验结果,并能掌握借助时域分析方法所求得的参数分析语音信号的基音周期及共振峰。

二、实验原理及实验结果
1.窗口的选择 通过对发声机理的认识,语音信号可以认为是短时平稳的。在 5~50ms 的范围内,语音频谱特性和一些物 理特性参数基本保持不变。我们将每个短时的语音称为一个分析帧。一般帧长取 10~30ms。我们采用一个长度 有限的窗函数来截取语音信号形成分析帧。 通常会采用矩形窗和汉明窗。 图 1.1 给出了这两种窗函数在帧长 N=50 时的时域波形。
矩 形 窗 2 1.8 1.6 1.4 1.2 1 0.9 0.8 0.7 0.6 hanming窗

w( n)

1 0.8 0.6 0.4 0.2 0

w( n)
0 20 sample 40 60

0.5 0.4 0.3 0.2 0.1 0

0

20 sample

40

60

图1.1

矩形窗和Hamming窗的时域波形

矩形窗的定义:一个N点的矩形窗函数定义为如下
?n? N w( n ) ? ? 1,0 0,其他

hamming窗的定义:一个N点的hamming窗函数定义为如下

? 0.54?0.46cos(2? Nn?1),0? n? N w(n)= ? 0,其他 ?
这两种窗函数都有低通特性,通过分析这两种窗的频率响应幅度特性可以发现(如图 1.2) :矩形窗的主瓣 宽度小(4*pi/N) ,具有较高的频率分辨率,旁瓣峰值大(-13.3dB) ,会导致泄漏现象;汉明窗的主瓣宽 8*pi/N, 旁瓣峰值低(-42.7dB) ,可以有效的克服泄漏现象,具有更平滑的低通特性。因此在语音频谱分析时常使用汉 明窗,在计算短时能量和平均幅度时通常用矩形窗。表 1.1 对比了这两种窗函数的主瓣宽度和旁瓣峰值。

矩形窗频率响应 0 -20

幅 度 /dB

-40 -60 -80

0

0.1

0.2

0.3

0.4 0.5 0.6 归 一 化 频 率 (f/fs) Hamming窗 频 率 响 应

0.7

0.8

0.9

1

0

幅 度 /dB

-50

-100

0

0.1

0.2

0.3

0.4 0.5 0.6 归 一 化 频 率 (f/fs)

0.7

0.8

0.9

1

图1.2

矩形窗和Hamming窗的频率响应

表1.1 矩形窗和hamming窗的主瓣宽度和旁瓣峰值

窗函数 矩形窗 hamming 2.短时能量

主瓣宽度 4*pi/N 8*pi/N

旁瓣峰值 13.3dB 42.7dB

由于语音信号的能量随时间变化,清音和浊音之间的能量差别相当显著。因此对语音的短时能量进行 分析,可以描述语音的这种特征变化情况。定义短时能量为:

En ?

m ???

? [ x(m)w(n ? m)]2 ?

?

m ? n ? N ?1

?

n

[ x(m)w(n ? m)]2
,其中 N 为窗长

特殊地,当采用矩形窗时,可简化为:

En ?

m ???

? x ( m)
2

?

图 1.3 和图 1.4 给出了不同矩形窗和 hamming 窗长的短时能量函数,我们发现:在用短时能量反映语音信 号的幅度变化时,不同的窗函数以及相应窗的长短均有影响。hamming 窗的效果比矩形窗略好。但是,窗的长 短影响起决定性作用。窗过大(N 很大) ,等效于很窄的低通滤波器,不能反映幅度 En 的变化;窗过小( N 很 小) ,短时能量随时间急剧变化,不能得到平滑的能量函数。在 11.025kHz 左右的采样频率下,N 选为 100~200 比较合适。 短时能量函数的应用:1)可用于区分清音段与浊音段。En 值大对应于浊音段,En 值小对应于清音段。2) 可用于区分浊音变为清音或清音变为浊音的时间(根据 En 值的变化趋势) 。3)对高信噪比的语音信号,也可以 用来区分有无语音(语音信号的开始点或终止点)。无信号(或仅有噪声能量)时,En 值很小,有语音信号时, 能量显著增大。

采样幅度

0 -1 0 2000 4000 6000 8000 sample 10000 12000 14000 16000 18000

采样幅度

1

1 0 -1 0 2000 4000 6000 8000 sample 10000 12000 14000 16000 18000

短时能量

2 0 0 2000 4000 6000 8000 sample 10000 12000 14000 16000

N=50 18000

短时能量

4

2 1 0 0 2000 4000 6000 8000 sample 10000 12000 14000 16000 N=50 18000

短时能量

5 0 0 2000 4000 6000 8000 sample 10000 12000 14000 16000

N=150 18000

短时能量

10

4 2 0 0 2000 4000 6000 8000 sample 10000 12000 14000 16000 N=150 18000

短时能量

5 0 0 2000 4000 6000 8000 sample 10000 12000 14000 16000

N=250 18000

短时能量

10

10 5 0 0 2000 4000 6000 8000 sample 10000 12000 14000 16000 N=250 18000

短时能量

10 0 0 2000 4000 6000 8000 sample 10000 12000 14000 16000

N=350 18000

短时能量

20

10 5 0 0 2000 4000 6000 8000 sample 10000 12000 14000 16000 N=350 18000

短时能量

10 0 0 2000 4000 6000 8000 sample 10000 12000 14000 16000

N=450 18000

短时能量

20

10 5 0 0 2000 4000 6000 8000 sample 10000 12000 14000 16000 N=450 18000

图 1.3 不同矩形窗长的短时能量函数

图 1.4

不同 hamming 窗长的短时能量函数

3.短时平均过零率 过零率可以反映信号的频谱特性。当离散时间信号相邻两个样点的正负号相异时,我们称之为“过零”,即 此时信号的时间波形穿过了零电平的横轴。统计单位时间内样点值改变符号的次数具可以得到平均过零率。定 义短时平均过零率:

Zn ?

m ???

?

?

sgn[ x[m] ? sgn[ x(m ? 1)] w(n ? m)

其中 sgn[] 为符号函数,

x ( n )?0 sgn x ( n) ? ? 1, ?1, x ( n ) 0

,在矩形窗条件下,可以简化为

Zn ?

1 2N

m ? n ? N ?1

?

n

sgn[ x(m) ? sgn[ x(m ? 1)]

短时过零率可以粗略估计语音的频谱特性。由语音的产生模型可知,发浊音时,声带振动,尽管声道有多 个共振峰,但由于声门波引起了频谱的高频衰落,因此浊音能量集中于 3KZ 以下。而清音由于声带不振动,声 道的某些部位阻塞气流产生类白噪声,多数能量集中在较高频率上。高频率对应着高过零率,低频率对应着低 过零率,那么过零率与语音的清浊音就存在着对应关系。. 图 1.5 为某一语音在矩形窗条件下求得的短时能量和短时平均过零率。分析可知:清音的短时能量较低, 过零率高,浊音的短时能量较高,过零率低。清音的过零率为 0.5 左右,浊音的过零率为 0.1 左右,但两者分布 之间有相互交叠的区域,所以单纯依赖于平均过零率来准确判断清浊音是不可能的,在实际应用中往往是采用 语音的多个特征参数进行综合判决。 短时平均过零率的应用:1)区别清音和浊音。例如,清音的过零率高,浊音的过零率低。此外,清音和浊 音的两种过零分布都与高斯分布曲线比较吻合。2)从背景噪声中找出语音信号。语音处理领域中的一个基本问 题是,如何将一串连续的语音信号进行适当的分割,以确定每个单词语音的信号,亦即找出每个单词的开始和 终止位置。3)在孤立词的语音识别中,可利用能量和过零作为有话无话的鉴别。

1

采样幅度

0.5

0

-0.5

0

2000

4000

6000

8000 sample

10000

12000

14000

16000

18000

8 6
短时能量

4 2 0

0

2000

4000

6000

8000 sample

10000

12000

14000

16000

18000

0.5
短时平均过零率

0.4 0.3 0.2 0.1 0 0 2000 4000 6000 8000 sample 10000 12000 14000 16000 18000

图 1.5 矩形窗条件下的短时平均过零率

4、短时自相关函数 自相关函数用于衡量信号自身时间波形的相似性。清音和浊音的发声机理不同,因而在波形上也存在着较 大的差异。浊音的时间波形呈现出一定的周期性,波形之间相似性较好;清音的时间波形呈现出随机噪声的特 性,样点间的相似性较差。因此,我们用短时自相关函数来测定语音的相似特性。短时自相关函数定义为:

Rn (k ) ?

m ???
'?

?

?

x(m) w(n ? m) x(m ? k ) w(n ? m ? k )
'

令 m ? n ? m ,并且 w(?m) ? w (m) ,可以得到:

Rn (k ) ?

m ???

? [ x(n ? m)w' (m)][ x(n ? m ? k )w' (m ? k )] ?

?

N ?1? k m?0

? [ x(n ? m)w (m)][ x(n ? m ? k )w (m ? k )]
' '

图 6 给出了清音的短时自相关函数波形,图 7 给出了不同矩形窗长条件下(窗长分别为 N=70,N=140,N=210, N=280)浊音的短时自相关函数波形。由图 1.6、图 1.7 短时自相关函数波形分析可知:清音接近于随机噪声, 清音的短时自相关函数不具有周期性,也没有明显突起的峰值,且随着延时 k 的增大迅速减小;浊音是周期信 号,浊音的短时自相关函数呈现明显的周期性,自相关函数的周期就是浊音信号的周期,根据这个性质可以判 断一个语音信号是清音还是浊音,还可以判断浊音的基音周期。浊音语音的周期可用自相关函数中第一个峰值 的位置来估算。所以在语音信号处理中,自相关函数常用来作以下两种语音信号特征的估计: 1)区分语音是清音还是浊音; 2)估计浊音语音信号的基音周期。
0.08 0.06 0.04 0.02 0 -0.02 -0.04 -0.06 -0.08 0 50 100 150 200 250 300

清音 0.1

0.05

R(k)

0

-0.05

-0.1

0

50

100

150 延时k

200

250

300

图 1.6

清音的短时自相关函数

5 N=70

R(k)

0

-5

0

20

40

60

80

100 延时k

120

140

160

180

200

220

5 N=140

R(k)

0

-5

0

20

40

60

80

100 延时k

120

140

160

180

200

220

10 N=210

R(k)

0

-10

0

20

40

60

80

100 延时k

120

140

160

180

200

220

10 N=280

R(k)

0

-10

0

20

40

60

80

100 延时k

120

140

160

180

200

220

图 1.7

不同矩形窗长条件下的浊音的短时自相关函数

5、时域分析方法的应用 1)基音频率的估计 首先可利用时域分析(短时能量、短时过零率、短时自相关)方法的某一个特征或某几个特征的结合,判 定某一语音有效的清音和浊音段;其次,针对浊音段,可直接利用短时自相关函数估计基音频率,其方法是: 估算浊音段第一最大峰的位置,再利用抽样率计算基音频率,举例来说,若某一语音浊音段的第一最大峰值约 为 35 个抽样点,设抽样频率为 11.025KHZ,则基音频率为 11025/35=315 HZ。 但是,实际上第一最大峰值位置有时并不一定与基音周期吻合。一方面与窗长有关,另一方面还与声道特 性有关。鉴于此,可采用三电平削波法先进行预处理。 2)语音端点的检测与估计 可利用时域分析(短时能量、短时过零率、短时自相关)方法的某一个特征或某几个特征的结合,判定某 一语音信号的端点,尤其在有噪声干扰时,如何准确检测语音信号的端点,这在语音处理中是富有挑战性的一 个课题。

三、附录(参考程序)
1) 短时能量 (1)加矩形窗 a=wavread('beifeng.wav'); subplot(6,1,1),plot(a); N=32; for i=2:6 h=linspace(1,1,2.^(i-2)*N);%形成一个矩形窗,长度为 2.^(i-2)*N En=conv(h,a.*a);% 求短时能量函数En subplot(6,1,i),plot(En); if(i==2) legend('N=32'); elseif(i==3) legend('N=64');

elseif(i==4) legend('N=128'); elseif(i==5) legend('N=256'); elseif(i==6) legend('N=512'); end end (2)加汉明窗 a=wavread('beifeng.wav'); subplot(6,1,1),plot(a); N=32; for i=2:6 h=hanning(2.^(i-2)*N);%形成一个汉明窗,长度为 2.^(i-2)*N En=conv(h,a.*a);% 求短时能量函数En subplot(6,1,i),plot(En); if(i==2) legend('N=32'); elseif(i==3) legend('N=64'); elseif(i==4) legend('N=128'); elseif(i==5) legend('N=256'); elseif(i==6) legend('N=512'); end end 2) 短时平均过零率 a=wavread('beifeng.wav'); n=length(a); N=320; subplot(3,1,1),plot(a); h=linspace(1,1,N); En=conv(h,a.*a); %求卷积得其短时能量函数 En subplot(3,1,2),plot(En); for i=1:n-1 if a(i)>=0 b(i)= 1; else b(i) = -1; end if a(i+1)>=0 b(i+1)=1; else b(i+1)= -1; end w(i)=abs(b(i+1)-b(i)); %求出每相邻两点符号的差值的绝对值 end k=1;

j=0; while (k+N-1)<n Zm(k)=0; for i=0:N-1; Zm(k)=Zm(k)+w(k+i); end j=j+1; k=k+N/2; %每次移动半个窗 end for w=1:j Q(w)=Zm(160*(w-1)+1)/(2*N); %短时平均过零率 end subplot(3,1,3),plot(Q),grid;

3) 自相关函数 N=240 Y=WAVREAD('beifeng.wav'); x=Y(13271:13510); x=x.*rectwin(240); R=zeros(1,240); for k=1:240 for n=1:240-k R(k)=R(k)+x(n)*x(n+k); end end j=1:240; plot(j,R); grid;

实验四
一、实验目的

基于 MATLAB 分析语音信号频域特征

信号的傅立叶表示在信号的分析与处理中起着重要的作用。因为对于线性系统来说,可以很方便地确定其 对正弦或复指数和的响应,所以傅立叶分析方法能完善地解决许多信号分析和处理问题。另外,傅立叶表示使 信号的某些特性变得更明显,因此,它能更深入地说明信号的各项物理现象。 由于语音信号是随着时间变化的,通常认为,语音是一个受准周期脉冲或随机噪声源激励的线性系统的输 出。输出频谱是声道系统频率响应与激励源频谱的乘积。声道系统的频率响应及激励源都是随时间变化的,因 此一般标准的傅立叶表示虽然适用于周期及平稳随机信号的表示,但不能直接用于语音信号。由于语音信号可 以认为在短时间内,近似不变,因而可以采用短时分析法。 本实验要求掌握傅里叶分析原理,会利用已学的知识,编写程序估计短时谱、倒谱,画出语谱图,并分析 实验结果,在此基础上,借助频域分析方法所求得的参数分析语音信号的基音周期或共振峰。

二、实验原理
1、短时傅立叶变换 由于语音信号是短时平稳的随机信号,某一语音信号帧的短时傅立叶变换的定义为:

X n (e jw ) ?

m???

? x(m)w(n ? m)e

?

? jwm

(2.1)

其中 w(n-m)是实窗口函数序列,n 表示某一语音信号帧。令 n-m=k',则得到

X n (e jw ) ?
于是可以得到

k '???

? w(k ') x(n ? k ')e
?

?

? jw( n ?k ')

(2.2)

X n (e jw ) ? e? jwn
假定

k ???

? w(k ) x(n ? k )e

jwk

(2.3)

X n (e jw ) ?
则可以得到

k ???

? w(k ) x(n ? k )e

?

jwk

(4)

X n (e jw ) ? e? jwn X n (e jw )

(5)

同样,不同的窗口函数,将得到不同的傅立叶变换式的结果。由上式可见,短时傅立叶变换有两个变量:n 和 ω,所以它既是时序 n 的离散函数,又是角频率 ω 的连续函数。与离散傅立叶变换逼近傅立叶变换一样,如 令 ω=2πk/N,则得离散的短时傅立叶吧如下:

Xn(e j 2? k / N ) ? Xn(k ) ?
m???

? x(m)w(n ? m)e

?

? j 2? km / N

,(0 ? k ? N ? 1)
(6)

2、语谱图 水平方向是时间轴,垂直方向是频率轴,图上的灰度条纹代表各个时刻的语音短时谱。语谱图反映了语音 信号的动态频率特性,在语音分析中具有重要的实用价值。被成为可视语言。 语谱图的时间分辨率和频率分辨率是由窗函数的特性决定的。时间分辨率高,可以看出时间波形的每个周 期及共振峰随时间的变化,但频率分辨率低,不足以分辨由于激励所形成的细微结构,称为宽带语谱图;而窄 带语谱图正好与之相反。 宽带语谱图可以获得较高的时间分辨率,反映频谱的快速时变过程;窄带语谱图可以获得较高的频率分辨 率,反映频谱的精细结构。两者相结合,可以提供带两与语音特性相关的信息。语谱图上因其不同的灰度,形

成不同的纹路,称之为“声纹”。声纹因人而异,因此可以在司法、安全等场合得到应用。

3、复倒谱和倒谱 复倒谱 x(n) 是 x(n)的 Z 变换取对数后的逆 Z 变换,其表达式如下:
^

x?Z

^

?1

[ln Z [ x(n)]]

(7)

倒谱 c(n)定义为 x(n)取 Z 变换后的幅度对数的逆 Z 变换,即

c(n) ? z ?1[ln | X ( z) |]

(8)

在时域上,语音产生模型实际上是一个激励信号与声道冲激响应的卷积。对于浊音,激励信号可以由周期 脉冲序列表示;对于清音,激励信号可以由随机噪声序列表示。声道系统相当于参数缓慢变化的零极点线性滤 波器。这样经过同态处理后,语音信号的复倒谱,激励信号的复倒谱,声道系统的复倒谱之间满足下面的关系:

s ( n ) ? e( n ) ? v ( n )

^

^

^

(9)

由于倒谱对应于复倒谱的偶部,因此倒谱与复倒谱具有同样的特点,很容易知道语音信号的倒谱,激励信 号的倒谱以及声道系统的倒谱之间满足下面关系:

c (n) ? c (n) ? c (n)
s e v

(10)

浊音信号的倒谱中存在着峰值,它的出现位置等于该语音段的基音周期,而清音的倒谱中则不存在峰值。 利用这个特点我们可以进行清浊音的判断,并且可以估计浊音的基音周期。

4、基因周期估计 浊音信号的倒谱中存在峰值,它的出现位置等于该语音段的基音周期,而清音的倒谱中则不存在峰值。利 用倒谱的这个特点,我们可以进行语音的清浊音判决,并且可以估计浊音的基音周期。首先计算语音的倒谱, 然后在可能出现的基因周期附近寻找峰值。如果倒谱峰值超过了预先设置的门限,则输入语音判断为浊音,其 峰值位置就是基因周期的估计值;反之,如果没有超出门限的峰值的话,则输入语音为清音。

5、共振峰估计 对倒谱进行滤波,取出低时间部分进行进行逆特征系统处理,可以得到一个平滑的对数谱函数,这个对数 谱函数显示了输入语音段的共振峰结构,同时谱的峰值对应于共振峰频率。通过此对数谱进行峰值检测,就可 以估计出前几个共振峰的频率和强度。对于浊音的声道特性,可以采用前三个共振峰来描述;清音不具备共振 峰特点。

三、实验结果
1 短时谱

original signal 1 0.5 0 -0.5 -1

0

2

4

6 短时谱

8

10

12 x 10
4

50

0

-50

-100

0

50

100

150

200

250

300

图 2.1 短时谱 2 语谱图

图 2.2

语谱图

3 倒谱和复倒谱 图 3、4 是加矩形窗和汉明窗的倒谱图和复倒谱图,图中横轴的单位是 Hz,纵轴的单位是 dB。

加矩形窗时的倒谱 1 0.5 0 -0.5 -1

0

50

100

150

200

250

300

加矩形窗时的复倒谱 5

0

-5

0

50

100

150

200

250

300

图 2.4 加矩形窗时的倒谱和复倒谱图

加汉明窗时的倒谱 1

0

-1

-2

0

50

100

150

200

250

300

加汉明窗时的复倒谱 20 10 0 -10 -20

0

50

100

150

200

250

300

图 2.3 加汉明窗时倒谱和复倒谱图

4 基因周期和共振峰估计

1 0

倒谱幅度

-1 -2 -3

0

100

200

300 点数N

400

500

600

100

幅 度 /dB

0

-100

-200

0

100

200

300 时 间 /ms

400

500

600

图 2.5 倒谱图 分析第 15 帧其中第一峰值出现在第 2 个样点,窗长为 512(64ms) ,抽样频率为 11KHz,说明基因频率就 在这个点上,其基因频率为 5.5KHz,基音周期为 0.182ms。 四、附录(参考程序) 1)短时谱
clear a=wavread('beifeng.wav'); subplot(2,1,1), plot(a);title('original signal'); grid N=256; h=hamming(N); for m=1:N b(m)=a(m)*h(m) end y=20*log(abs(fft(b))) subplot(2,1,2) plot(y);title('短时谱'); grid

2)语谱图 [x,fs,nbits]=wavread('beifeng.wav')
specgram(x,512,fs,100); xlabel('时间(s)'); ylabel('频率(Hz)');

title('语谱图');

3)倒谱和复倒谱 (1)加矩形窗时的倒谱和复倒谱
clear a=wavread('beifeng.wav',[4000,4350]); N=300; h=linspace(1,1,N); for m=1:N b(m)=a(m)*h(m); end c=cceps(b); c=fftshift(c); d=rceps(b); d=fftshift(d); subplot(2,1,1) plot(d);title('加矩形窗时的倒谱') subplot(2,1,2) plot(c);title('加矩形窗时的复倒谱')

(2)加汉明窗时的倒谱和复倒谱
clear a=wavread('beifeng.wav',[4000,4350]); N=300; h=hamming(N); for m=1:N b(m)=a(m)*h(m); end c=cceps(b); c=fftshift(c); d=rceps(b); d=fftshift(d); subplot(2,1,1) plot(d);title('加汉明窗时的倒谱') subplot(2,1,2) plot(c);title('加汉明窗时的复倒谱')

实验五
一、实验目的

基于 MATLAB 的 LPC 分析

线性预测分析是最有效的语音分析技术之一,在语音编码、语音合成、语音识别和说话人识别等语音处理 领域中得到了广泛的应用。语音线性预测的基本思想是:一个语音信号的抽样值可以用过去若干个取样值的线 性组合来逼近。通过使实际语音抽样值与线性预测抽样值的均方误差达到最小,可以确定唯一的一组线性预测 系数。 采用线性预测分析不仅能够得到语音信号的预测波形,而且能够提供一个非常好的声道模型。如果将语音 模型看作激励源通过一个线性时不变系统产生的输出,那么可以利用 LP 分析对声道参数进行估值,以少量低 信息率的时变参数精确地描述语音波形及其频谱的性质。此外,LP 分析还能够对共振峰、功率谱等语音参数进 行精确估计,LP 分析得到的参数可以作为语音识别的重要参数之一。 由于语音是一种短时平稳信号,因此只能利用一段语音来估计模型参数。此时有两种方案:一种是将长的 语音序列加窗,然后对加窗语音进行 LP 分析,只要限定窗的长度就可以保证分析的短时性,这种方案称为自 相关法;另一种方案不对语音加窗,而是在计算均方预测误差时限制其取和区间,这样可以导出 LP 分析的自 协方差法。 本实验要求掌握 LPC 原理,会利用已学的知识,编写程序估计线性预测系数以及 LPC 的推演参数,并能 利用所求的相关参数估计语音的端点、清浊音判断、基因周期、共振峰等。

二、实验原理
1 LP 分析基本原理 LP 分析为线性时不变因果稳定系统 V(z)建立一个全极点模型,并利用均方误差准则,对已知的语音信 号 s(n)进行模型参数估计。 如果利用 P 个取样值来进行预测, 则称为 P 阶线性预测。 假设用过去 P 个取样值 的加权之和来预测信号当前取样值
S ? n ? ? ? ak ? n ? k ?
k ?1 ? p

?S ? n ?1? , S ? n ? 2? ,

S ? n ? p ??

S ? n?

,则预测信号

S ?n?

?

为:

(1)

其中加权系数用 a k 表示,称为预测系数,则预测误差为:
e ? n ? ? s ? n ? ? S ? n ? ? s ? n ? ? ? ak ? n ? k ?
k ?1 ? p

(2)

要使预测最佳,则要使短时平均预测误差最小有:
2 ? ? E? ?e ? n ? ? ? ? min
2 ?? ?e ? n ? ? ?

(3)

?ak

? 0, (1 ? k ? p)

(4)



? ?i, k ? ? E ? ? s ? n ? i ? , S ? n ? k ?? ?
最小的 ? 可表示成:

(5)

? min ? ? ? 0,0? ? ? ak? ? 0, k ?
k ?1

p

(6)

显然,误差越接近于零,线性预测的准确度在均方误差最小的意义上为最佳,由此可以计算出预测系数。 通过 LPC 分析,由若干帧语音可以得到若干组 LPC 参数,每组参数形成一个描绘该帧语音特征的矢量, 即 LPC 特征矢量。由 LPC 特征矢量可以进一步得到很多种派生特征矢量,例如线性预测倒谱系数、线谱对特 征、部分相关系数、对数面积比等等。不同的特征矢量具有不同的特点,它们在语音编码和识别领域有着不同 的应用价值。 2 自相关法 在最佳线性预测中,若用下式定义的时间平均最小均方准则代替(3)式的集合平均最小均方准则,即令

??

1 N

N ? p ?1 n ?0

?

e2 ? n ? ? min

(7)

事实上就是短时自相关函数,因而
R ? i ? k ? ? ? ? i, k ?

(8) (9)

R ?k ? ? E ? ? S ? n ? , S ? n ? k ?? ?

根据平稳随机信号的自相关性质,可得

? ?i, k ? ? R ? i ? k ? , i ? 1, 2
由(6)式,可得:

p; k ? 0,1

p

(10)

? min ? R ? 0? ? ? ak R ? k ?
k ?1

p

(11)

综上所述,可以得到如下矩阵形式:
R ?1? ? R ? 0? ? R ?0? ? R ?1? ? ? ? R ? P ? 1? R ? P ? 2 ? ? ? ? R ? P ? 1? ? ? R ? P ? 2? ? ? ? R ?0? ? ? ? ?

? a1 ? ? R ?1? ? ? ? ? R ? 2? ? ? ? a2 ? ? ? a3 ? ? ? R ? 3? ? ? ? ? ? ? ? ? ? ? ?a ? ? R p ? ? n? ? ? ? ??

(12)

值得注意的是,自相关法在计算预测误差时,数据段

?S ?0? , S ?1? ,

S ? n ?1??

的两端都需要加 P 个零取样值,

因而可造成谱估计失真。特别是在短数据段的情况下,这一现实更为严重。另外,当预测系数量化时,有可能 造成实际系统的不稳定。

自相关解法主要有杜宾算法、格型算法和舒尔算法等几种高效递推算法。 3 协方差法 如果在最佳线性预测中,用下式定义的时间平均最小均方准则代替(3)式的集合平均最小均方准则,则可得 到类似的方程:
??
1 N
N ?1 n? p

? e2 ? n ? ? min

(13)

可以看出,这里的数据段两端不需要添加零取样值。在理论上,协方差法计算出来的预测系数有可能造成 预测误差滤波器的不稳定,但在实际上当每帧信号取样足够多时,其计算结果将与自相关法的结果很接近,因 而稳定性一般是能够保证的 (当然这种方法也有量化效应可能引起不稳定的缺点)。 协方差解法的最大优点在于不存在自相关法中两端出现很大预测误差的情况,在 N 和 P 相差不大时,其参 数估值比自相关法要精确的多。但是在语音信号处理时,往往取 N 在 200 左右。此时,自相关法具有较大误差 的段落在整个语音段中所占的比例很小,参数估值也是比较准确的。在这种情况下,协方差法误差较小的优点 就不再突出,其缺乏高效递推算法的缺点成为了制约因素。所以,在语音信号处理中往往使用高效的自相关法。 4 全极点声道模型 将线性预测分析应用于语音信号处理,不仅是为了利用其预测功能,更因为它提供了一个非常好的声道模 型。 将式(2)所示的方程看成是滤波器在语音信号激励下的输入输出方程,则该滤波器称为预测误差滤波器,其 e(n)是输出误差。变换到 z 域,P 阶预测误差滤波器的系统函数为
i H ? z? ? 1 ? ? ia ?z i ?1 p

(14)

可以看出,如果将预测误差 e(n)作为激励信号,使其通过预测误差滤波器的逆滤波器 H(Z),即
H ? z? ? 1 ? A? Z ? 1 1 ? ? ai z ?i
i ?1 p

(15)
ai ? i ? 1, 2, , p?

则 H(Z)的输出为语音信号 s(n),也就是说,H(Z)在预测误差 e(n)的激励下可以合成语音。因此,H(Z)被称 为语音信号的全极点模型,也称为语音合成器。该模型的参数就是 P 阶线性预测的预测系数 。

因为预测误差含有语音信号的基音信息,所以对于浊音,模型的激励信号源是以基音周期重复的单位脉冲; 对于清音,激励信号源 e(n)是自噪声。语音信号的全极点模型是一种很重要的声道模型,是许多应用和研究的 基础。 5 LPCC 如果声道特性 H(Z)用式(14)所示的全极点模型表示,有
H ? z? ? S ? z? I ? z? ? 1 1 ? ? an z ? n
n ?1 p

(16)

式中,S(z)和 I(z)分别为语音信号 s n 和激励源 in 的 Z 变换。对人的听觉来说,浊音是最重要的语音信号。对于浊 音,模型的激励信号源 e(n)是以基音周期重复的单位脉冲,此时有
I ? z? ? 1

。可得 s n 的 Z 变换 S(z)为
S ?z? ? 1 1 ? ? an z ? n
n ?1 p

(17)

式中,

ai ? i ? 1, 2,

, p?

为 P 阶线性预测系数。根据倒谱的定义,对具有最小相位特征的语音信号 s n ,有
?

ln S ? z ? ? C ? z ? ? ? cn z ?n
n ?1

(18)
?1

式中, c n 为语音信号的倒谱。将式(16)代入式(17),并对两边 z 求导,得
?c1 ? a1 ? n ?1 ? ? k? ?cn ? an ? ? ?1 ? n ? ak cn ? k ,1 ? n ? p ? k ?1 ? ?

(19)

根据上式即可由线性预测系数通过递推得到倒谱系数,将这样得到的倒谱称为线性预测倒谱系数。 6 结合语音帧能量构成LPC组合参数 由于人能从声音的音色、频高等各种信息中感知说话人的个性,因此可以想象,利用特征的有效组合可以 得到比较稳定的识别性能。一般来说,如果组合的各参量之间相关性不大,则会更有效一些,因为它们分别反 映了语音信号中的不同特征。多年来,人们对组合参数在说话人识别中的应用进行了大量研究 。实验证明,组 合参数可以提高系统的识别性能。 组合参数虽然可以提高系统的性能,但很显然,无论是在特征参数提取环节,还是在模型训练和模型匹配 环节都使运算量有所增加。在特征参数提取环节,要计算一种以上的特征参数。在模型训练和模型匹配环节, 由于组合参数特征矢量的维数较多,使运算复杂度有所增加。运算量的增加会使系统的识别速度受到影响。 为使运算量问题得到较好的解决,所以可以由 LPC 参数与语音帧能量构成组合参数,能够在运算量增加不 明显的情况下改进系统的性能。 语音帧能量是指一帧语音信号的能量,它等于该帧语音样值的平方和。选取与语音帧能量构成组合参数主 要有以下考虑:1)语音帧能量是语音信号最基本的短时参数之一,它表征一帧语音信号能量的大小,是语音信 号一个重要的时域特征;2)由一帧语音求出的语音帧能量是一个标量值,与其它参量构成组合参数不会使原特 征矢量的维数明显增加,特征矢量的维数越少,则需要的运算复杂度越小,另外,获取语音帧能量的运算并不 复杂;3)语音帧能量与 LPC 参数之间的相关性不大,它们反映的是语音信号的不同特征,应该有较好的效果。 7 模型增益G 模型的激励信号
p

Ge ? n ?

表示为: (20)

Ge ? n ? ? s ? n ? ? ? ai s ? n ? i ?
i ?1

预测误差 e(n)如式(2) ,这样当实际的预测系数与模型系数相等时,有

? ? n ? ? Ge ? n ?

(21)

这说明激励信号正比于误差信号,其比例常数等于模型增益 G。通常假设误差信号的能量等于输入激励信 号的能量,因此可以得到:
G2 ? e2 ? m? ? ? ? 2 ? m? ? En
m?0 m?0 N ?1 N ?1

(22)

对于式中的激励信号

e ? n?

,主要分为浊音和清音两种情况。其中为浊音时,考虑到此时实际的激励信号为声门

脉冲,因此可以将激励信号表示为 n ? 0 时的单位抽样。为了保证这个假设成立,要求分析的区间应该大致和语 音基因周期的长度相等。当语音为清音时,我们假定激励信号 采用自相关解法时,浊音的模型增益为
En ? Rn ? 0? ? ? ai Rn ?i ? ? G2
i ?1 p

e ? n?

为一个零均值、单位方差的平稳白噪声过程。

(23)

清音计算模型增益的公式和浊音相同。

三、实验结果(参考)
我们使用的原始语音为“北风”,采样频率为 11000Hz,运行程序见附录。 在这里我们取第 30 帧进行观察,线性预测阶数为 12,看到图 3.1 所示的原始语音帧的波形,预测语音帧波 形和它们之间预测误差的波形。图 3.2 为原始语音帧和预测语音帧的短时谱和 LPC 谱的波形
原始语音波形 1 0 -1

0

2

4

6

8

10

12 x 10
4

原始语音和预测语音波形 0.5 0 -0.5

0

50

100

150 预测误差

200

250

300

0.2 0 -0.2

0

50

100

150

200

250

300

图3.1 原始语音帧、预测语音帧和预测误差的波形

短时谱 100

0

幅度
-100 -200

0

10

20

30 40 频 率 /dB LPC谱

50

60

70

200 150

幅度

100 50 0

0

10

20

30 40 频 率 /dB

50

60

70

图3.2 原始语音帧和预测语音帧的短时谱和LPC谱的波形

这里我们可以改变线性误差的阶数来观察语音帧的短时谱和LP谱的变化情况,如图3.3。
P1 =5 100

0

幅度
-100 -200

0

10

20

30 频 率 /dB P1 =10

40

50

60

70

100

0

幅度
-100 -200

0

10

20

30 频 率 /dB P1 =20

40

50

60

70

100

0

幅度
-100 -200

0

10

20

30 频 率 /dB

40

50

60

70

图3.3 预测阶数对语音帧短时谱和LPC谱的影响

从图中可以看出,P 越大,LPC 谱越能反映出语音短时谱的细节部分,但 LPC 谱的光滑度随之下降。由于 我们的目的只是用 LPC 谱反映声道综合效应的谱的表示式,而具体的谐波形状是通过激励谱来控制的,因此 LPC 谱只要能够体现出语音的共振峰的结构和谱包络就可以,因此从计算复杂性的角度分析,预测阶数 P 应该 适中。 图 3.4 是原始语音和预测误差的倒谱波形,我们可以从中计算出原始语音的基音周期。从图中看出两峰值 之间的间隔为 40 点左右,基音周期为 40/11000=3.6ms,频率为 278Hz 左右。

原始语音帧倒谱 1

0

/dB
-1 -2 0

50

100

150 语音帧 预测误差倒谱

200

250

300

1

0

/dB
-1 -2 0

50

100

150 语音帧

200

250

300

图3.4 原始语音和预测误差的倒谱波形

图 3.5 给出了原始语音的语谱图和预测语音的语谱图,通过比较发现,预测语音的预测效果还可以,基音 频率相差无几。
原始语音语谱图 60

Frequency

40 20 0

0

100

200

300

400 500 600 Time 预测语音语谱图

700

800

900

60

Frequency

40 20 0

0

100

200

300

400 500 Time

600

700

800

900

图3.5 原始语音的语谱图和预测语音的语谱图

三、附录(LPC 分析参考程序) MusicSource = wavread('bei'); Music_source = MusicSource'; N = 256; % window length,N = 100 -- 1000; Hamm = hamming(N); % create Hamming window frame = input('请键入想要处理的帧位置 = '); % origin is current frame origin = Music_source(((frame - 1) * (N / 2) + 1):((frame - 1) * (N / 2) + N)); Frame = origin .* Hamm'; % %Short Time Fourier Transform %

[s1,f1,t1] = specgram(MusicSource,N,N/2,N); [Xs1,Ys1] = size(s1); for i = 1:Xs1 FTframe1(i) = s1(i,frame); end N1 = input('请键入预测器阶数 = '); % N1 is predictor's order [coef,gain] = lpc(Frame,N1); % LPC analysis using Levinson-Durbin recursion est_Frame = filter([0 -coef(2:end)],1,Frame); % estimate frame(LP) FFT_est = fft(est_Frame); err = Frame - est_Frame; % error % FFT_err = fft(err); subplot(2,1,1),plot(1:N,Frame,1:N,est_Frame,'-r');grid;title('原始语音帧vs.预测后语音帧') subplot(2,1,2),plot(err);grid;title('误差'); pause %subplot(2,1,2),plot(f',20*log(abs(FTframe2)));grid;title('短时谱') % % Gain solution using G^2 = Rn(0) - sum(ai*Rn(i)),i = 1,2,...,P % fLength(1 : 2 * N) = [origin,zeros(1,N)]; Xm = fft(fLength,2 * N); X = Xm .* conj(Xm); Y = fft(X , 2 * N); Rk = Y(1 : N); PART = sum(coef(2 : N1 + 1) .* Rk(1 : N1)); G = sqrt(sum(Frame.^2) - PART); A = (FTframe1 - FFT_est(1 : length(f1'))) ./ FTframe1 ; % inverse filter A(Z) subplot(2,1,1),plot(f1',20*log(abs(FTframe1)),f1',(20*log(abs(1 ./ A))),'-r');grid;title('短时谱'); subplot(2,1,2),plot(f1',(20*log(abs(G ./ A))));grid;title('LPC谱'); pause %plot(abs(ifft(FTframe1 ./ (G ./ A))));grid;title('excited') %plot(f1',20*log(abs(FFT_est(1 : length(f1')) .* A / G )));grid; %pause % % find_pitch % temp = FTframe1 - FFT_est(1 : length(f1')); % not move higher frequnce pitch1 = log(abs(temp)); pLength = length(pitch1);

result1 = ifft(pitch1,N); % move higher frequnce pitch1((pLength - 32) : pLength) = 0; result2 = ifft(pitch1,N); % direct do real cepstrum with err pitch = fftshift(rceps(err)); origin_pitch = fftshift(rceps(Frame)); subplot(211),plot(origin_pitch);grid;title('原始语音帧倒谱(直接调用函数)'); subplot(212),plot(pitch);grid;title('预测误差倒谱(直接调用函数)'); pause subplot(211),plot(1:length(result1),fftshift(real(result1)));grid;title('预测误差倒谱(根据定义编写,没有去除高频分 量)'); subplot(212),plot(1:length(result2),fftshift(real(result2)));grid;title('预测误差倒谱(根据定义编写,去除高频分量)');

实验六
一、实验目的

基于 VQ 的特定人孤立词语音识别研究

矢量量化(Vector Quantization)是一种极其重要的信号压缩方法,是自 70 年代末才发展起来的。它广泛应 用于语音编码、语音识别与合成、图象压缩等领域。VQ 在语音信号处理中占有十分重要的地位。许多重要的 研究课题中,特别是低速语音编码和语音识别的研究中,VQ 都起着非常重要的作用。 量化可以分为两大类:一类是标量量化,另一类是矢量量化。标量量化是将取样后的信号值逐个地进行量 化,而矢量量化是将若干个取样信号分成一组,即构成一个矢量,然后对此矢量一次进行量化。当然,矢量量 化压缩数据的同时也有信息的损失,但这仅取决于量化的精度。矢量量化是标量量化的发展,可以说,凡是要 用量化的地方都可以应用矢量量化。 本实验要求掌握矢量量化的原理,会利用已学的相关语音特征,构建语音特征矢量,然后利用 VQ 技术, 编写训练 VQ 码表的程序,并在此基础上利用所学的语音识别技术,编程实现基于矢量量化的特定人孤立词语 音识别,要注意的是识别过程中语音端点如何检测,从识别的实时性角度出发,建议能利用 VC 技术实现。

二、实验原理
1 矢量量化 1)基本原理 矢量量化的过程是:将语音信号波形的 K 个样点的每一帧,或者有 K 个参数的每一参数帧,构成 K 维空 间中的一个矢量,然后对这个矢量进行量化。通常所说的标量量化,也可以说是 K=1 的一维矢量量化。矢量量 化的过程与标量量化相似。在标量量化时,在一维的零至无大值之间设置若干个量化阶梯,当某输入信号的幅 度值落在某相邻的两个量化阶梯之间时,就被量化为两阶梯的中心值。而在矢量量化时,将 K 维无限空间划为 M 个区域边界,然后将输入矢量与这些边界进行比较,并被量化为“距离”最小的区域边界的中心矢量值。 2) 、失真测度 设计矢量量化器的关键是编码器的设计,而译码器的工作仅是一个简单的查表过程。在编码的过程中,需 要引入失真测度的概念。失真是将输入信号矢量用码书的重构矢量来表征时的误差或所付出的代价。而这种代 价的统计平均值(平均失真)描述了矢量量化器的工作特性。 在矢量量化器的设计中,失真测度的选择是很重要的。失真测度选用的合适与否,直接影响系统的性能。 要使所选择的失真测度有实际意义,必须具备以下几个条件:在主观评价上有意义,即最小的失真应该对应与 好的主观语言质量;易于处理,即在数学上易于实现,这样可以用于实际的矢量量化器的设计;平均失真存在 并且可以计算。 2 LBG 算法 算法是由 Linde,Buzo 和 Gray 在 1980 年首次提出的,常称为 LBG 算法。它是标量量化器中 Lioyd 算法的 多维推广。整个算法实际上就是反复迭代的过程,既用初始码书寻找最佳码书的迭代过程。它由对初始码书进 行迭代优化开始,一直到系统性能满足要求或者不再有明显的改进为止。 这种算法既可以用于已知信号源概率分布的场合,也可以用于未知信号源概率分布的场合,但此时要知道 它的一系列输出值(称为训练序列) 。由于通常语音信号的概率分布随着各种应用场合的不同,不可能事先统计

过,因而无法知道它的概率分布。所以目前多用训练序列来设计码书和矢量量化器。 3 语音识别 语音识别是研究使机器能够准确地听出人的语音内容的问题,即准确的识别所说的语音。语音识别是近二 三十几年发展起来的新兴学科,在计算机、信息处理、通信与电子系统、自动控制等领域中有着广泛的应用。 运用语音识别技术,人们设计了各种语音识别系统。有的已经应用于实际,有的还处在研究阶段。其中对 孤立词的识别,研究的最早也最成熟,目前,对孤立词的识别无论是小词汇量还是大词汇量,无论是与讲话者 有关还是与讲话者无关,在实验室中的正识率已经达到 95%以上。 这种系统存在的问题最少,因为单词之间有停顿,可以使识别问题简单化;且单词之间的端点检测比较容 易;单词之间的协同发音影响也可以减至最低;对孤立词的发音都比较认真。由于此系统本身用途广泛,且其 许多技术对其他类型系统有通用性并易于推广,所以稍加补充一些知识就可用于其他类型系统(如在识别部分 加用适当语义信息等,则可用于连续语音识别中) 。采用矢量量化技术主要用于减少计算量,应用于特征处理可 减少特征的类型从而减少计算量,也可以推广应用到摸板的归并压缩。其主要工作就是聚类,即在特征空间中 合理的拟定一组点(称为一组聚类中心或码本) ,每个中心称为码字。于是特征空间中任一点均可按最小距离准 则用码本之一来代表。 不管用何种语音识别方法,主要过程由两部分组成,一是训练,一是识别。在进行训练时,用观察的序列 训练得到参考模型集,每一个模型对应于摸板中的一个单词。在进行识别时,为每一个参考模型计算出产生测 试观察的概率,且测试信号(即输入信号)按最大概率被识别为某个单词。 要实现上面的隐马尔可夫模型,模型的输入信号必须取自有限字母集中的离散序列,也就是说,必须将连 续的语音信号变为有限离散的序列。例如若模型的输入信号为 LPC 参数这样的矢量信号,那么用矢量量化完成 上述的识别过程是非常合适的。下面是 VQ/HMM 孤立单词识别的方框图。图中矢量量化器作为整个识别系统 的一个前处理器。

三、实验结果
就算法模型方面而言,需要有进一步的突破。目前能看出它的一些明显不足,尤其在中文语音识别方面,语言 模型还有待完善,因为语言模型和声学模型正是听写识别的基础,这方面没有突破,语音识别的进展就只能是一句 空话。 目前使用的语言模型只是一种概率模型,还没有用到以语言学为基础的文法模型,而要使计算机确实理解人 类的语言,就必须在这一点上取得进展,这是一个相当艰苦的工作。此外,随着硬件资源的不断发展,一些核心算法 如特征提取、搜索算法或者自适应算法将有可能进一步改进。可以相信,半导体和软件技术的共同进步将为语音 识别技术的基础性工作带来福音。 就强健性方面而言,语音识别技术需要能排除各种环境因素的影响。 目前,对语音识别效果影响最大的就是环 境杂音或嗓音,在公共场合,你几乎不可能指望计算机能听懂你的话,来自四面八方的声音让它茫然而不知所措。 很显然这极大地限制了语音技术的应用范围 , 目前 , 要在嘈杂环境中使用语音识别技术必须有特殊的抗嗓 (NoiseCancellation)麦克风才能进行,这对多数用户来说是不现实的。 在公共场合中,个人能有意识地摒弃环境嗓音 并从中获取自己所需要的特定声音,如何让语音识别技术也能达成这一点呢?这的确是一个艰巨的任务。 以下是实验部分的截图:

程序附录
using System; using System.Collections.Generic; using System.ComponentModel; using System.Data; using System.Drawing; using System.Text; using System.Windows.Forms; using SpeechLib; namespace WindowsApplication1 {public partial class Form1 : Form {private SpVoice voice; //语音合成对象 private SpSharedRecoContext objRecoContext; //共享型上下文对象 private ISpeechRecoGrammar grammar; //上下文语法对象添加成员函数 public void RecoContext_Recongnition(int StreamNumber, object StreamPosition, SpeechRecognitionType RecongitionType, ISpeechRecoResult Result) { textBox1.Text += Result.PhraseInfo.GetText(0, -1, true) + ""; //输出语音识别结果 } private void InitializeSpeech() {//语音识别初始化 objRecoContext = new SpSharedRecoContext();//初始化共享型上下文对象 objRecoContext.Recognition += new _ISpeechRecoContextEvents_RecognitionEventHandler(RecoContext_Recongnition);//为上下文添加事件处理函数 grammar = objRecoContext.CreateGrammar(0);//创建语法器 grammar.DictationLoad("", SpeechLoadOption.SLOStatic);//装载语法规则 grammar.DictationSetState(SpeechLib.SpeechRuleState.SGDSActive);//激活语法器 objRecoContext.State = SpeechRecoContextState.SRCS_Disabled; }

public Form1() { InitializeComponent(); } private void button1_Click(object sender, EventArgs e) {objRecoContext.State = SpeechRecoContextState.SRCS_Disabled; voice.Speak(textBox1.Text, SpeechVoiceSpeakFlags.SVSFlagsAsync); } private void Form1_Load(object sender, EventArgs e) { voice = new SpVoiceClass(); this.InitializeSpeech(); } private void button2_Click(object sender, EventArgs e) {objRecoContext.State = SpeechRecoContextState.SRCS_Enabled; } private void button3_Click(object sender, EventArgs e) { objRecoContext.State = SpeechRecoContextState.SRCS_Disabled; } } } 参考文献 1: 基于 VQ 的特定人孤立词语音识别, 徐海儿 参考文献 2: 基于 VQ 的语音识别系统, 谢福香 参考文献 3: 基于 VQ 的特定人孤立词语音识别,高虹 参考文献 4:基于 VQ 的特定人孤立词识别,徐志浩 参考文献 5:基于 VQ 的特定人孤立词语音识别,来海静

实验七
1、实验目的 (1)了解回声的产生和梳妆滤波器; (2)混音效果的原理和均衡器的设计; 2、 实验要求

音乐信号处理

(1) 设计函数实现一段语音或音乐的回声产生; (2)设计均衡器,使得得不同频率的混合音频信号,通过一个均衡器后,增强或削减某些频率区 域,以便修正低频和高频信号之间的关系; 3、 实验提示 ( 1 )回声产生可以使用梳妆滤波器, y(n)=x(n)+ax(n-R), a<1( 回声衰减系数 ) ;或者传输函数为
H ( z) ?

? ? z ?R ,? ?1 1 ? ?z ? R 的全通滤波器实现;比较这两种实现方式的区别,分析为什么会有这

样的区别; (2) 可以用许多一阶和二阶参数可调的滤波器级联来实现均衡器的功能,滤波器的结构选择结构 要求是调整方便,最好调一个参数只影响一个应用指标,且可调参数少;

实验八
1、实验目的

双音多频(DTMF)信号的合成和识别

了解电话按键音形成的原理,理解 DTMF 音频产生软件和 DTMF 解码算法; 利用 FFT 算法识别按键音; 2、 实验要求 设计音频产生函数,音频信号见下图,每个数据信号持续半秒; 实现解码函数:接受(1)产生的 DTMF 信号,识别信号的频率,并生成包含拨号数字的 序列; 3、 实验提示 DTFT 音频可以用两个正弦波按比例叠加产生; 检测算法可以用 FFT 算法的 DFT,或是用一组滤波器实现; Goertzel 算法可以实现调谐滤波器;

附页 1: 1 函数 mffc: function r = mfcc(s, fs) % s 声音信号的向量 fs 取样频率 m = 100; n = 256; l = length(s); nbFrame = floor((l - n) / m) + 1; for i = 1:n for j = 1:nbFrame M(i, j) = s(((j - 1) * m) + i); end end h = hamming(n); M2 = diag(h) * M; for i = 1:nbFrame frame(:,i) = fft(M2(:, i)); end t = n / 2; tmax = l / fs; m = melfb(20, n, fs); n2 = 1 + floor(n / 2); z = m * abs(frame(1:n2,: )).^2; r = dct(log(z)); 2 主程序: 17 / 18 % Demo script that generates all graphics in the report and demonstrates our results. [s1 fs1] = wavread('s1.wav'); [s2 fs2] = wavread('s2.wav'); %Question 2 disp('> Question 2:画出原始语音波形'); t = 0:1/fs1:(length(s1) - 1)/fs1; plot(t, s1), axis([0, (length(s1) - 1)/fs1 -0.4 0.5]); title('原始语音 s1 的波形'); xlabel('时间/s'); ylabel('幅度') pause close all %Question 3 (linear) disp('> Question 3: 画出线性谱'); M = 100;%当前帧数 N = 256;%帧长 frames = blockFrames(s1, fs1, M, N);%分帧 t = N / 2;

tm = length(s1) / fs1; subplot(121); imagesc([0 tm], [0 fs1/2], abs(frames(1:t, :)).^2), axis xy; title('能量谱(M = 100, N = 256)'); xlabel('时间/s'); ylabel('频率/Hz'); colorbar; %Question 3 (logarithmic) disp('> Question 3: 画出对数谱'); subplot(122); imagesc([0 tm], [0 fs1/2], 20 * log10(abs(frames(1:t, :)).^2)), axis xy; title('对数能量谱(M = 100, N = 256)'); xlabel('时间/s'); ylabel('频率/Hz'); colorbar; D=get(gcf,'Position'); set(gcf,'Position',round([D(1)*.5 D(2)*.5 D(3)*2 D(4)*1.3])) pause close all %Question 4 disp('> Question 4: 画出不同帧长语谱图'); lN = [128 256 512]; u=220; for i = 1:length(lN) N = lN(i); M = round(N / 3); frames = blockFrames(s1, fs1, M, N); t = N / 2; temp = size(frames); nbframes = temp(2); u=u+1; subplot(u) imagesc([0 tm], [0 fs1/2], 20 * log10(abs(frames(1:t, :)).^2)), axis xy; title(sprintf(' 能量 对数谱 ( 第 = %i 帧 , 帧长 = %i, 帧数 = %i)', M, N, nbframes)); xlabel('时间/s'); ylabel('频率/Hz'); colorbar end D=get(gcf,'Position'); set(gcf,'Position',round([D(1)*.5 D(2)*.5 D(3)*1.5 D(4)*1.5])) pause close all %Question 5

disp('> Question 5: Mel 空间'); plot(linspace(0, (fs1/2), 129), (melfb(20, 256, fs1))'); title('Mel 滤波'); xlabel('频率/Hz'); pause close all %Question 6 disp('> Question 6: 修正谱'); M = 100; N = 256; frames = blockFrames(s1, fs1, M, N); n2 = 1 + floor(N / 2); m = melfb(20, N, fs1); z = m * abs(frames(1:n2, :)).^2; t = N / 2; tm = length(s1) / fs1; subplot(121) imagesc([0 tm], [0 fs1/2], abs(frames(1:n2, :)).^2), axis xy; title('原始能量谱'); xlabel('时间/s'); ylabel('频率/Hz'); colorbar; subplot(122) imagesc([0 tm], [0 20], z), axis xy; title('通过 mel 倒谱修正后的能量谱'); xlabel('时间/s'); ylabel('滤波器数目'); colorbar; D=get(gcf,'Position'); set(gcf,'Position',[0 D(2) D(3)*2 D(4)]) pause close all %Question 7 disp('> Question 7: 2D plot of accustic vectors'); c1 = mfcc(s1, fs1); c2 = mfcc(s2, fs2); plot(c1(5, :), c1(6, :), 'or'); hold on; plot(c2(5, :), c2(6, :), 'xb'); xlabel('5th Dimension'); ylabel('6th Dimension'); legend('说话人 1', '说话人 2'); title('2D plot of accoustic vectors'); pause

close all %Question 8 disp('> Question 8: 画出已训练好的 VQ 码本') d1 = vqlbg(c1,16); d2 = vqlbg(c2,16); plot(c1(5, :), c1(6, :), 'xr') hold on plot(d1(5, :), d1(6, :), 'vk') plot(c2(5, :), c2(6, :), 'xb') plot(d2(5, :), d2(6, :), '+k') xlabel('5th Dimension'); ylabel('6th Dimension'); legend('说话人 1', '码本 1', '说话人 2', '码本 2'); title('2D plot of accoustic vectors'); pause close all 3 函数 blockFrames function M3 = blockFrames(s, fs, m, n) % blockFrames: % Puts the signal into frames 分帧函数 % Inputs: % s contains the signal to analize 语音信号 % fs is the sampling rate of the signal 语音采样频率 % m is the distance between the beginnings of two frames 两帧之间的距 离 % n is the number of samples per frame 每帧的采样点数 % Output: % M3 is a matrix containing all the frames 数组形式,包含了所有的帧 l = length(s); %语音信号的长度 nbFrame = floor((l - n) / m) + 1; %帧数 for i = 1:n for j = 1:nbFrame M(i, j) = s(((j - 1) * m) + i); %逐帧扫描 end end h = hamming(n); M2 = diag(h) * M; %加汉明窗 for i = 1:nbFrame M3(:, i) = fft(M2(:, i)); %短时傅立叶变换 End 4 函数 disteu function d = disteu(x, y)

% DISTEU Pairwise Euclidean distances between columns of two matrices 测试失 真度 % Input: % x, y: Two matrices whose each column is an a vector data. % Output: % d: Element d(i,j) will be the Euclidean distance between two column vectors X(:,i) and Y(:,j) % Note: % The Euclidean distance D between two vectors X and Y is: % D = sum((x-y).^2).^0.5 [M, N] = size(x); [M2, P] = size(y); if (M ~= M2) error('不匹配!') end d = zeros(N, P); if (N < P) copies = zeros(1,P); for n = 1:N d(n,:) = sum((x(:, n+copies) - y) .^2, 1); end else copies = zeros(1,N); for p = 1:P d(:,p) = sum((x - y(:, p+copies)) .^2, 1)'; end end d = d.^0.5; 5.函数 vqlbg function r = vqlbg(d,k) % VQLBG Vector quantization using the Linde-Buzo-Gray algorithme 使用 LBG 算法 % % Inputs: % d contains training data vectors (one per column) % k is number of centroids required % Output: % r contains the result VQ codebook (k columns, one for each centroids) e = .01; r = mean(d, 2); dpr = 10000; for i = 1:log2(k) r = [r*(1+e), r*(1-e)];

while (1 == 1) z = disteu(d, r); [m,ind] = min(z, [], 2); t = 0; for j = 1:2^i r(:, j) = mean(d(:, find(ind == j)), 2); x = disteu(d(:, find(ind == j)), r(:, j)); for q = 1:length(x) t = t + x(q); end end if (((dpr - t)/t) < e) break; else dpr = t; end end end 6 函数 test function test(testdir, n, code) % Speaker Recognition: Testing Stage % Input: % testdir : string name of directory contains all test sound files % n : number of test files in testdir 音频文件数目 % code : codebooks of all trained speakers 说话人的训练码本 % Note: % Sound files in testdir is supposed to be: % s1.wav, s2.wav, ..., sn.wav % Example: % >> test('C:\data\test\', 8, code); for k = 1:n % 读出音频文件 file = sprintf('%ss%d.wav', testdir, k); [s, fs] = wavread(file); v = mfcc(s, fs); distmin = inf; k1 = 0; for l = 1:length(code) % 计算每个训练码本的失真度 d = disteu(v, code{l}); dist = sum(min(d,[],2)) / size(d,1); % 计算 MFCC's

if dist < distmin distmin = dist; k1 = l; end end msg = sprintf('Speaker %d matches with speaker %d', k, k1); disp(msg); end 7 函数 train function code = train(traindir, n) % 计算 wav 文件的 VQ 码码本 % Speaker Recognition: Training Stage % Input: % traindir : string name of directory contains all train sound files % n : number of train files in traindir % Output: % code : trained VQ codebooks, code{i} for i-th speaker % Note: % Sound files in traindir is supposed to be: % s1.wav, s2.wav, ..., sn.wav % Example: % >> code = train('C:\data\train\', 8); k = 16; % number of centroids required for i = 1:n % train a VQ codebook for each speaker file = sprintf('%ss%d.wav', traindir, i); disp(file); [s, fs] = wavread(file); v = mfcc(s, fs); % Compute MFCC's code{i} = vqlbg(v, k); % Train VQ codebook end 8、函数 melfb function m = melfb(p, n, fs) % MELFB Determine matrix for a mel-spaced filterbank % Inputs: p number of filters in filterbank 滤波器数 % n length of fft FFT 变换的点数 % fs sample rate in Hz 采样频率 % Outputs: x a (sparse) matrix containing the filterbank amplitudes % size(x) = [p, 1+floor(n/2)] % Usage: For example, to compute the mel-scale spectrum of a % colum-vector signal s, with length n and sample rate fs: %

% f = fft(s); % m = melfb(p, n, fs); % n2 = 1 + floor(n/2); % z = m * abs(f(1:n2)).^2; % % z would contain p samples of the desired mel-scale spectrum % % To plot filterbanks e.g.: % % plot(linspace(0, (12500/2), 129), melfb(20, 256, 12500)'), % title('Mel-spaced filterbank'), xlabel('Frequency (Hz)'); f0 = 700 / fs; fn2 = floor(n/2); lr = log(1 + 0.5/f0) / (p+1); % convert to fft bin numbers with 0 for DC term bl = n * (f0 * (exp([0 1 p p+1] * lr) - 1)); b1 = floor(bl(1)) + 1; b2 = ceil(bl(2)); b3 = floor(bl(3)); b4 = min(fn2, ceil(bl(4))) - 1; pf = log(1 + (b1:b4)/n/f0) / lr; fp = floor(pf); pm = pf - fp; r = [fp(b2:b4) 1+fp(1:b3)]; c = [b2:b4 1:b3] + 1; v = 2 * [1-pm(b2:b4) pm(1:b3)]; m = sparse(r, c, v, p, 1+fn2)
附录2:
1.分帧函数:
function M3 = blockFrames(s, fs, m, n) l = length(s); %语音信号的长度 nbFrame = floor((l - n) / m) + 1; for i = 1:n for j = 1:nbFrame M(i, j) = s(((j - 1) * m) + i); end end h = hamming(n); M2 = diag(h) * M; for i = 1:nbFrame %加汉明窗 %逐帧扫描 %帧数

M3(:, i) = fft(M2(:, i)); End 2.失真测度函数: function d = disteu(x, y) [M, N] = size(x); [M2, P] = size(y); if (M ~= M2) error('不匹配!') end d = zeros(N, P); if (N < P) copies = zeros(1,P); for n = 1:N

%短时傅立叶变换

d(n,:) = sum((x(:, n+copies) - y) .^2, 1); end else copies = zeros(1,N); for p = 1:P d(:,p) = sum((x - y(:, p+copies)) .^2, 1)'; end end d = d.^0.5; 3.三角滤波器组 function m = melfb(p, n, fs) f0 = 700 / fs; fn2 = floor(n/2); lr = log(1 + 0.5/f0) / (p+1); bl = n * (f0 * (exp([0 1 p p+1] * lr) - 1)); b1 = floor(bl(1)) + 1; b2 = ceil(bl(2)); b3 = floor(bl(3)); b4 = min(fn2, ceil(bl(4))) - 1; pf = log(1 + (b1:b4)/n/f0) / lr; fp = floor(pf); pm = pf - fp; r = [fp(b2:b4) 1+fp(1:b3)]; c = [b2:b4 1:b3] + 1; v = 2 * [1-pm(b2:b4) pm(1:b3)]; m = sparse(r, c, v, p, 1+fn2); 4.MFCC function r = mfcc(s, fs) m = 100; n = 256; l = length(s); nbFrame = floor((l - n) / m) + 1; for i = 1:n%分帧 % s声音信号的向量,fs取样频率

for j = 1:nbFrame M(i, j) = s(((j - 1) * m) + i); end end h = hamming(n);%加窗 M2 = diag(h) * M; for i = 1:nbFrame frame(:,i) = fft(M2(:, i));%离散傅立叶变换 end t = n / 2; tmax = l / fs; m = melfb(20, n, fs);%通过一组Mel尺度的三角形滤波器组 n2 = 1 + floor(n / 2); z = m * abs(frame(1:n2, :)).^2;%频谱幅度的平方,得到能量谱 r = dct(log(z));%计算每个滤波器组输出的对数能量,再经过经过离散余弦反变换,得到MFCC系数 5.LBG算法 function r = vqlbg(d,k) e = .01;%阈值 r = mean(d, 2);%取提取出来的所有帧的特征矢量的型心(均值)作为第一个码字矢量(初始码书) dpr = 10000;%算法最大迭代次数 for i = 1:log2(k)%分裂算法 r = [r*(1+e), r*(1-e)]; while (1 == 1) z = disteu(d, r);%训练矢量量化失真量的总和及相对失真 [m,ind] = min(z, [], 2); t = 0; for j = 1:2^i%重新计算各个区域的新型心,得到新的码书 r(:, j) = mean(d(:, find(ind == j)), 2); x = disteu(d(:, find(ind == j)), r(:, j)); for q = 1:length(x) t = t + x(q); end end if (((dpr - t)/t) < e)%若相对失真小于阈值 break;%迭代结束 else dpr = t; end end end 6.训练码本 function code = train(traindir, n) k = 16; % 质心数目 % 训练码本 for i = 1:n

file = sprintf('%ss%d.wav', traindir, i); disp(file);

[s, fs] = wavread(file); v = mfcc(s, fs); end 7.测试语音 function test(testdir, n, code) for k = 1:n % 读出音频文件 file = sprintf('%ss%d.wav', testdir, k); [s, fs] = wavread(file); v = mfcc(s, fs); % 计算MFCC's distmin = inf; k1 = 0; for l = 1:length(code) d = disteu(v, code{l}); dist = sum(min(d,[],2)) / size(d,1); % 计算说话人的平均量化失真 if dist < distmin distmin = dist; k1 = l; end end msg = sprintf('说话人 %d 和说话人 %d语音相匹配', k, k1); disp(msg); end 8.总体演示(demo) [s1 fs1] = wavread('s1.wav'); [s2 fs2] = wavread('s2.wav'); %Question 2 disp('> Question 2:画出原始语音波形'); t = 0:1/fs1:(length(s1) - 1)/fs1; plot(t, s1), axis([0, (length(s1) - 1)/fs1 -0.4 0.5]); title('原始语音s1的波形'); xlabel('时间/s'); ylabel('幅度') pause close all %Question 3 (linear) disp('> Question 3: 画出线性谱'); M = 100;%当前帧数 N = 256;%帧长 frames = blockFrames(s1, fs1, M, N);%分帧 t = N / 2; tm = length(s1) / fs1; % 计算MFCC code{i} = vqlbg(v, k); % 生成的VQ码本

subplot(121); imagesc([0 tm], [0 fs1/2], abs(frames(1:t, :)).^2), axis xy; title('能量谱(M = 100, N = 256)'); xlabel('时间/s'); ylabel('频率/Hz'); colorbar; %Question 3 (logarithmic) disp('> Question 3: 画出对数谱'); subplot(122); imagesc([0 tm], [0 fs1/2], 20 * log10(abs(frames(1:t, :)).^2)), axis xy; title('对数能量谱(M = 100, N = 256)'); xlabel('时间/s'); ylabel('频率/Hz'); colorbar; D=get(gcf,'Position'); set(gcf,'Position',round([D(1)*.5 D(2)*.5 D(3)*2 D(4)*1.3])) pause close all %Question 4 disp('> Question 4: 画出不同帧长语谱图'); lN = [128 256 512]; u=220; for i = 1:length(lN) N = lN(i); M = round(N / 3); frames = blockFrames(s1, fs1, M, N); t = N / 2; temp = size(frames); nbframes = temp(2); u=u+1; subplot(u) imagesc([0 tm], [0 fs1/2], 20 * log10(abs(frames(1:t, :)).^2)), axis xy; title(sprintf('能量对数谱(第 = %i帧, 帧长 = %i, 帧数 = %i)', M, N, nbframes)); xlabel('时间/s'); ylabel('频率/Hz'); colorbar end D=get(gcf,'Position'); set(gcf,'Position',round([D(1)*.5 D(2)*.5 D(3)*1.5 D(4)*1.5])) pause close all %Question 5 disp('> Question 5: Mel空间'); plot(linspace(0, (fs1/2), 129), (melfb(20, 256, fs1))'); title('Mel滤波');

xlabel('频率/Hz'); pause close all %Question 6 disp('> Question 6: 修正谱'); M = 100; N = 256; frames = blockFrames(s1, fs1, M, N); n2 = 1 + floor(N / 2); m = melfb(20, N, fs1); z = m * abs(frames(1:n2, :)).^2; t = N / 2; tm = length(s1) / fs1; subplot(121) imagesc([0 tm], [0 fs1/2], abs(frames(1:n2, :)).^2), axis xy; title('原始能量谱'); xlabel('时间/s'); ylabel('频率/Hz'); colorbar; subplot(122) imagesc([0 tm], [0 20], z), axis xy; title('通过mel倒谱修正后的能量谱'); xlabel('时间/s'); ylabel('滤波器数目'); colorbar; D=get(gcf,'Position'); set(gcf,'Position',[0 D(2) D(3)*2 D(4)]) pause close all %Question 7 disp('> Question 7: 二维矢量图'); c1 = mfcc(s1, fs1); c2 = mfcc(s2, fs2); plot(c1(5, :), c1(6, :), 'or'); hold on; plot(c2(5, :), c2(6, :), 'xb'); xlabel('5th Dimension'); ylabel('6th Dimension'); legend('说话人1', '说话人2'); title('二维矢量图'); pause close all %Question 8 disp('> Question 8: 画出已训练好的VQ码本') d1 = vqlbg(c1,16); d2 = vqlbg(c2,16);

plot(c1(5, :), c1(6, :), 'xr') hold on plot(d1(5, :), d1(6, :), 'vk') plot(c2(5, :), c2(6, :), 'xb') plot(d2(5, :), d2(6, :), '+k') xlabel('5th Dimension'); ylabel('6th Dimension'); legend('说话人1', '码本1', '说话人2', '码本2'); title('二维矢量图'); pause close all


赞助商链接
更多相关文档:

数字语音信号处理实验

语音信号处理实验 2015 年 10 月 28 日 语音信号处理实验实验学时数:8 实验...通过开放实验, 目的使学生进一步理解数字语音信息处理的基本方法, 提高学生自主分...

语音信号处理实验指导书_图文

语音信号处理实验指导书 - 数字语音信号处理 实验指导书 编写 曹建荣 山东建筑大学信息与电气工程学院 2011 年 10 月 1 前言 语音信号处理是研究用数字信号处理...

数字语音信号处理实验__图文

数字语音信号处理 实验指导书 前言 语音信号处理是研究用数字信号处理技术和语音学知识对语音信号进行处理的新兴的 学科, 是目前发展最为迅速的信息科学研究领域的...

数字语音信号处理实验

数字语音信号处理实验_IT/计算机_专业资料。数字语音信号处理实验学号:1228401060 姓名:唐榆 专业:信息工程 日期:2014 年 12 月 21 日 目录 前言--- 数字语音信...

数字语音信号处理

数字语音信号处理 实验指导书 编著 王让定 熊益群 徐国娟 宁波大学信息科学与工程学院 2008 年 6 月 前言 语音信号处理是研究用数字信号处理技术和语音学知识对...

数字信号处理实验指导

数字语音信号处理 实验指导书 编著 王让定 熊益群 徐国娟 宁波大学信息科学与工程学院 2008 年 6 月 前言 语音信号处理是研究用数字信号处理技术和语音学知识对...

数字语音信号处理课程设计

( 论文 ) 题 目 数字语音信号处理 的程序设计 学生姓名 学号 林琳 0706020137...的数字语音信号分析原理说明 4.1.matlab 环境简介 MATLAB 是矩阵实验室(Matrix ...

数字语音信号处理

(设计)数字语音信号处理 学生姓名 指导教师 院、系、中心 专业年级 论文答辩...3 数字语音信号处理 而这项指标在 20 世纪 90 年代中后期实验室研究中得到了...

语音信号处理 实验 数字信号处理 傅里叶变换 DTFT

语音信号处理 实验 数字信号处理 傅里叶变换 DTFT_工学_高等教育_教育专区。语音信号处理 实验 数字信号处理 傅里叶变换 DTFT实验二 语音信号处理 语音信号处理综合...

《语音信号处理》期末试题总结

P62 9、语音信号处理也可以简称为语音处理,它是利用数字信号处理技术对语音信号进行处理 的一门学科,包括语音编码、 语音合成 、 语音识别 、说话人识别和 语音...

更多相关标签:
网站地图

文档资料共享网 nexoncn.com copyright ©right 2010-2020。
文档资料共享网内容来自网络,如有侵犯请联系客服。email:zhit325@126.com