数字语音信号处理实验（学生） - 图文

数字语音信号处理

实验指导书

北方学院信息科学与工程学院

电子教研室

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时的时域波形。

矩形窗21.81.61.41.210.90.80.70.6hanming窗w（n）10.80.60.40.20w（n）020sample40600.50.40.30.20.10020sample4060图1.1 矩形窗和Hamming窗的时域波形

矩形窗的定义：一个N点的矩形窗函数定义为如下

?n?Nw(n)??1,00,其他

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)

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、短时傅立叶变换

由于语音信号是短时平稳的随机信号，某一语音信号帧的短时傅立叶变换的定义为：

Xn(e)?

jwm????x(m)w(n?m)e??jwm （2.1）

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

于是可以得到

Xn(e)?jwk'????w(k')x(n?k')e???jw(n?k') （2.2）

Xn(e)?e 假定

jw?jwnk????w(k)x(n?k)ejwk （2.3）

Xn(e)? 则可以得到

jwk????w(k)x(n?k)e?jwk （4）

jw?jwnjwX(e)?eX(e) （5） nn

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

Xn(ej2?k/N)?Xn(k)?

2、语谱图

水平方向是时间轴，垂直方向是频率轴，图上的灰度条纹代表各个时刻的语音短时谱。语谱图反映了语音信号的动态频率特性，在语音分析中具有重要的实用价值。被成为可视语言。

语谱图的时间分辨率和频率分辨率是由窗函数的特性决定的。时间分辨率高，可以看出时间波形的每个周期及共振峰随时间的变化，但频率分辨率低，不足以分辨由于激励所形成的细微结构，称为宽带语谱图；而窄带语谱图正好与之相反。

宽带语谱图可以获得较高的时间分辨率，反映频谱的快速时变过程；窄带语谱图可以获得较高的频率分辨率，反映频谱的精细结构。两者相结合，可以提供带两与语音特性相关的信息。语谱图上因其不同的灰度，形

m????x(m)w(n?m)e??j2?km/N,(0?k?N?1) (6)

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

3、复倒谱和倒谱

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

^x?Z^?1[lnZ[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) （10）

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

4、基因周期估计

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

5、共振峰估计

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

三、实验结果

1 短时谱

original signal10.50-0.5-10246短时谱500-50-10081012x 104050100150200250300

图2.1 短时谱 2 语谱图

图2.2 语谱图

3 倒谱和复倒谱

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

加矩形窗时的倒谱10.50-0.5-1050100150200250300加矩形窗时的复倒谱50-5050100150200250300

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

加汉明窗时的倒谱10-1-2050100150200250300加汉明窗时的复倒谱20100-10-20050100150200250300

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

4 基因周期和共振峰估计

10倒谱幅度-1-2-30100200300点数N400500600100幅度/dB0-100-2000100200300时间/ms400500600

图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)

其中加权系数用ak表示，称为预测系数，则预测误差为:

e?n??s?n??S?n??s?n???ak?n?k?k?1?p (2)

要使预测最佳，则要使短时平均预测误差最小有：

2??E?e??n????min (3)

2???e?n????ak?0,(1?k?p) (4)

令

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

(5)

?min???0,0???ak??0,k?k?1p (6)

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

在最佳线性预测中，若用下式定义的时间平均最小均方准则代替(3)式的集合平均最小均方准则，即令 1??NN?p?1n?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,1p (10)

?min?R?0???akR?k?k?1p (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???Rp?n??????

(12)

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

?S?0?,S?1?,S?n?1??的两端都需要加P个零取样值，

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

自相关解法主要有杜宾算法、格型算法和舒尔算法等几种高效递推算法。 3 协方差法

如果在最佳线性预测中，用下式定义的时间平均最小均方准则代替(3)式的集合平均最小均方准则，则可得到类似的方程：

1??NN?1n?p?e2?n??min （13）

可以看出，这里的数据段两端不需要添加零取样值。在理论上，协方差法计算出来的预测系数有可能造成预测误差滤波器的不稳定，但在实际上当每帧信号取样足够多时，其计算结果将与自相关法的结果很接近，因而稳定性一般是能够保证的 (当然这种方法也有量化效应可能引起不稳定的缺点)。

协方差解法的最大优点在于不存在自相关法中两端出现很大预测误差的情况，在N和P相差不大时，其参数估值比自相关法要精确的多。但是在语音信号处理时，往往取N在200左右。此时，自相关法具有较大误差的段落在整个语音段中所占的比例很小，参数估值也是比较准确的。在这种情况下，协方差法误差较小的优点就不再突出，其缺乏高效递推算法的缺点成为了制约因素。所以，在语音信号处理中往往使用高效的自相关法。 4 全极点声道模型

将线性预测分析应用于语音信号处理，不仅是为了利用其预测功能，更因为它提供了一个非常好的声道模型。

将式(2)所示的方程看成是滤波器在语音信号激励下的输入输出方程，则该滤波器称为预测误差滤波器，其e(n)是输出误差。变换到z域，P阶预测误差滤波器的系统函数为

iH?z??1??ia?zi?1p （14）

可以看出，如果将预测误差e(n)作为激励信号，使其通过预测误差滤波器的逆滤波器H(Z)，即

H?z??1?A?Z?11??aiz?ii?1p （15）

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

ai?i?1,2,,p?。

因为预测误差含有语音信号的基音信息，所以对于浊音，模型的激励信号源是以基音周期重复的单位脉冲；对于清音，激励信号源e(n)是自噪声。语音信号的全极点模型是一种很重要的声道模型，是许多应用和研究的基础。 5 LPCC

如果声道特性H(Z)用式(14)所示的全极点模型表示，有

H?z??S?z?I?z??11??anz?nn?1p （16）

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

ai?i?1,2, （17）

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

?lnS?z??C?z???cnz?nn?1 （18）

?1式中，cn为语音信号的倒谱。将式(16)代入式(17)，并对两边z求导，得

?c1?a1?n?1??k?c?a??n?1??akcn?k,1?n?p?nn?k?1?? （19）

根据上式即可由线性预测系数通过递推得到倒谱系数，将这样得到的倒谱称为线性预测倒谱系数。 6 结合语音帧能量构成LPC组合参数

由于人能从声音的音色、频高等各种信息中感知说话人的个性，因此可以想象，利用特征的有效组合可以得到比较稳定的识别性能。一般来说，如果组合的各参量之间相关性不大，则会更有效一些，因为它们分别反映了语音信号中的不同特征。多年来，人们对组合参数在说话人识别中的应用进行了大量研究 。实验证明，组合参数可以提高系统的识别性能。

组合参数虽然可以提高系统的性能，但很显然，无论是在特征参数提取环节，还是在模型训练和模型匹配环节都使运算量有所增加。在特征参数提取环节，要计算一种以上的特征参数。在模型训练和模型匹配环节，由于组合参数特征矢量的维数较多，使运算复杂度有所增加。运算量的增加会使系统的识别速度受到影响。

为使运算量问题得到较好的解决，所以可以由LPC参数与语音帧能量构成组合参数，能够在运算量增加不明显的情况下改进系统的性能。

语音帧能量是指一帧语音信号的能量，它等于该帧语音样值的平方和。选取与语音帧能量构成组合参数主要有以下考虑：1）语音帧能量是语音信号最基本的短时参数之一，它表征一帧语音信号能量的大小，是语音信号一个重要的时域特征；2）由一帧语音求出的语音帧能量是一个标量值，与其它参量构成组合参数不会使原特征矢量的维数明显增加，特征矢量的维数越少，则需要的运算复杂度越小，另外，获取语音帧能量的运算并不复杂；3）语音帧能量与LPC参数之间的相关性不大，它们反映的是语音信号的不同特征，应该有较好的效果。 7 模型增益G

模型的激励信号

pGe?n?表示为：

（20）

Ge?n??s?n???ais?n?i?i?1预测误差e(n)如式（2），这样当实际的预测系数与模型系数相等时，有

??n??Ge?n? （21）

这说明激励信号正比于误差信号，其比例常数等于模型增益G。通常假设误差信号的能量等于输入激励信号的能量，因此可以得到：

G2?e2?m????2?m??Enm?0m?0N?1N?1 （22）

对于式中的激励信号

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

脉冲，因此可以将激励信号表示为n?0时的单位抽样。为了保证这个假设成立，要求分析的区间应该大致和语音基因周期的长度相等。当语音为清音时，我们假定激励信号

采用自相关解法时，浊音的模型增益为

En?Rn?0???aiRn?i??G2i?1pe?n?为一个零均值、单位方差的平稳白噪声过程。

(23)

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

三、实验结果(参考)

我们使用的原始语音为“北风”，采样频率为11000Hz,运行程序见附录。

在这里我们取第30帧进行观察，线性预测阶数为12，看到图3.1所示的原始语音帧的波形，预测语音帧波形和它们之间预测误差的波形。图3.2为原始语音帧和预测语音帧的短时谱和LPC谱的波形

原始语音波形10-1024681012x 104原始语音和预测语音波形0.50-0.50.20-0.2050100150预测误差200250300050100150200250300

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

短时谱1000幅度-100-200010203040频率/dBLPC谱506070200150幅度100500010203040频率/dB506070

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

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

P1 =51000幅度-100-2000102030频率/dBP1 =10405060701000幅度-100-2000102030频率/dBP1 =20405060701000幅度-100-2000102030频率/dB40506070

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

从图中可以看出，P越大，LPC谱越能反映出语音短时谱的细节部分，但LPC谱的光滑度随之下降。由于我们的目的只是用LPC谱反映声道综合效应的谱的表示式，而具体的谐波形状是通过激励谱来控制的，因此LPC谱只要能够体现出语音的共振峰的结构和谱包络就可以，因此从计算复杂性的角度分析，预测阶数P应该适中。

图3.4是原始语音和预测误差的倒谱波形，我们可以从中计算出原始语音的基音周期。从图中看出两峰值之间的间隔为40点左右，基音周期为40/11000=3.6ms，频率为278Hz左右。

原始语音帧倒谱10/dB-1-2050100150语音帧预测误差倒谱20025030010/dB-1-2050100150语音帧200250300

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

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

原始语音语谱图60Frequency402000100200300400500600Time预测语音语谱图70080090060Frequency402000100200300400500Time600700800900

图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(\装载语法规则

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(回声衰减系数)；或者传输函数为

??z?RH(z)?,??11??z?R的全通滤波器实现；比较这两种实现方式的区别，分析为什么会有这样的区别；

（2）可以用许多一阶和二阶参数可调的滤波器级联来实现均衡器的功能，滤波器的结构选择结构

要求是调整方便，最好调一个参数只影响一个应用指标，且可调参数少；

实验八 双音多频（DTMF）信号的合成和识别

1、实验目的

了解电话按键音形成的原理，理解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); % 计算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('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) % 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);%通过一组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); % 计算MFCC code{i} = vqlbg(v, k); % 生成的VQ码本 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;

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滤波');

本文来源：https://www.bwwdw.com/article/iezr.html