南邮通达语音时频参数的提取和应用和语音编码实验.doc

1、南京邮电大学通达学院实验报告实验名称语音时频参数的提取和应用实验语音编码实验课程名称信息处理技术专业综合实验班级学号姓名开课时间2016/2017学年 ,第 二 学期实验一语音时频参数的提取和应用一、实验目的熟悉基本的语音时频参数提取方法，观察并比较各类参数采用不同提取方法的效果，了解这些时频参数在语音处理中的应用。二、实验内容1、编写并调试语音时域参数提取程序，包括短时功率、短时过零率等函数。2、编写并调试语音频域参数提取程序，包括频谱、语谱图、基音频率、共振峰等。3、利用提取出的时频参数对语音信号进行静/ 清 /浊音的划分。三、实验原理语音信号具有短时平稳性，需对其进行分段（帧）处理（帧长

2、一般取10 至 30ms之间）。根据激励方式的不同，语音信号可以分成浊音和清音两大类：浊音：激励源可模拟成准周期性的脉冲串。该周期称为基音周期，其倒数就是基音频率，基音频率一般在几百赫兹以内。清音：激励源可模拟成随机白噪声。简化的语音生成数学模型如图所示：基音周期准周期脉冲序列发生器浊音时变线性清音系统语音取样值随机噪声发生器增益控制时变参数对语音时 频波形进行绘制与观察，有助于了解语音的清期等声门参数和共振峰等声道参数进行提取与分析，/浊音等特性。对基音周是语音处理与应用的重要内容。对于时变语音信号x( n)，可以计算它的短时功率Px (m) 和短时过零率Z x (m) 。第m帧（帧长为N个

3、样点）的功率计算公式为：1m2x(n)Px ( m)Nn m N 1过零率计算公式为：1msign x(n) sign x(n 1)Z x ( m)2N n m N 1对于浊音语音，可以利用其频谱X () 具有丰富的谐波分量的特点，求出其谐波乘积谱，计算公式为：RHPSx()X ( r)r1式中， R一般取为5。在谐波乘积谱中，基频分量变得很大，更易于估计基音周期。四、实验方法及程序1. 调用 matlab 中的 wavread 和 wavplay 读入并播放语音文件。2. 调用 fft 命令计算语音的频谱。3. 调用 specgram 命令得到语音的语谱图。4. 调用 plot 命令绘制并观

4、察语音的时频域波形。5. 根据短时功率、短时过零率、谐波乘积谱的计算公式，编写分别用以计算短时 功 率 、 短 时 过 零 率 、 基 音 频 率 的 函 数 ： functionPx=stpower、 functionstzerocross(x, N)和 function HPSx = hpspectrum(x,N,R)。Zx=6. 利用短时功率和短时过零率这两个参数可以对语音信号进行浊音编写实现此分类功能的函数voiunvoi.m。/清音分类。五、实验结果与分析1. 清音和浊音的短时功率、短时过零率各有何特点？清音段的能量一般比浊音段的小得多。浊音时能量集中于较低频率段内，具有较低的过零率

5、，而清音时能量集中于较高频率段内，具有较高的过零率。2. 如何利用短时功率和短时过零率这两个参数对语音信号进行浊音/ 清音的分类？发浊音时，由于声门波引起了谱的高频跌落。所以其语音能量集中在3kHz以下；发清音时，多数能量出现在高频。高频率意味着高的平均过零率，低频率意味着低的平均过零率，一般可认为浊音具有较低的平均过零率，清音具有较高的平均过零率。3. 如何根据时域波形估计各帧元音的基音周期？看图可知四帧语音的周期均为0.01s ，第一帧0100Hz 有 10个峰值点，基频约为10Hz ，第二帧0100Hz基频约为20Hz ，第四帧有 5 个峰值点，基频约为20Hz ，第三帧0100Hz01

6、00Hz 有 4个峰值点，基频约为25Hz 。有 5 个峰值点，4. 如何根据基于 DFT 的对数幅度谱估计某幀浊音的基音周期？01172 为一个基因周期。5. 如何根据基于 DFT 的对数幅度谱，估计出共振峰频率？有 7个峰点，因此基频约为167 。6. 时域对语音信号进行加窗，反映在频域，其窗谱对基于DFT的对数幅度谱有何影响？如何估计出窗谱的主瓣宽度？因为时 间窗 幅度 的傅 立叶 变换 为脉 冲状 态 ,则频谱偏差 很小 , 当 定常 过程 x(n) 的平均值 为零且 N 很大时ETN()=f()(12)所以加窗的周期图TN( ) 是频谱密度函数的无偏估计。固有TN()=1+ ( )f

7、()。如果原始信号的频谱成份与FFT 中的谱线完全一致，这种情况下采样数据的长度为信号周期的整数倍，频谱中只有主瓣。7. 如何根据窄带、宽带语谱图，提取出语音的基音和共振峰轨迹？在基音提取中，广泛采用语音波形或误差信号波形的低通滤波。提取共振峰特性最简便的手段是使用语谱仪。提取共振峰还有倒谱法、析法等更为有效、准确的方法。（做完试验后，记住要回答这些问题）LPC分实验过程记录与结果分析1、进入 matlab1.1 在 Windows 环境中，创建一个名为speech 的文件夹，存放语音数据和与实验相关的 Matlab 文件。如：E:speech 。用于实验的语音数据（.mat ）及相应说明（.

8、txt ）包括：digits.mat%digits.txtgliss.mat%两个包含 /i/ 的滑音gliss.txtletters.mat%英文字母表中26 个字母的发音letters.txtma1.txtma1-1.mat%一句连续语音ma1-2.mat%一句连续语音ma1-3.mat%一句连续语音ma1-4.mat%一句连续语音timit.txttimit1.mat%一句连续语音timit2.mat%一句连续语音timit3.mat%一句连续语音timit4.mat%一句连续语音vowels.mat%元音 /a/, /i/, /o/, /u/ 的发音vowels.txtwords.ma

9、t%十个孤立字的发音words.txt1.2在 Windows桌面上，双击Matlab的图标，进入Matlab令窗口 提示符后，键入cd E:speech指令，进入实验文件夹路径。的工作环境。在命1.3 在命令窗口键入clear all指令，清除工作空间的所有变量。2、加载语音数据，熟悉matlab 命令2.1 加载“ timit1.mat”语音数据，命令为：load timit12.2用命令“who ”和“whos ”列出当前工作空间中的变量。变量“timit1”会以包含55911 个元素的矩阵（列向量）形式出现Your variables are:timit1NameSizeBytesCl

10、assAttributestimit155911x1447288double2.3 利用命令“length”和“size ”查看语音信号timit1的长度和维数：m = length(timit1)m,n = size(timit1)（粘贴实验结果）m =55911m =55911n =12.4 通过以下命令，分别查看语音信号timit1中间的1 个、 2 个和 20 个数据值：timit1(5001)timit1(5001,5003)timit1(5001:5020)（粘贴实验结果）ans =-0.0646ans =-0.06460.0111ans =-0.06460.03070.0111-

11、0.02000.20210.40120.1869-0.1756-0.19020.08230.23630.17120.0637-0.0209-0.0683-0.03940.00810.01050.03110.06962.5 Matlab中，符号“”可以实现转置运算，运行以下命令进行观察：size (timit1)（粘贴实验结果）ans =1559112.6 使用函数max 和 min 可以分别得到信号的最大值和最小值：max(timit1)min(timit1)（粘贴实验结果）ans =1ans =-0.80683、语音信号的时域分析3.1 使用以下命令，画出语音信号timit1 （采样频率为1

12、6kHz ，长度为3.5 秒）的时域波形。plot (timit1)title ( timit1 )xlabel ( Time Index,it n)ylabel ( Amplitude)axis (0 length(timit1) -1 1)grid3.2 分别画出一帧浊音和一帧清音的语音时域波形（采样频率为25ms ，每帧有400 个样点） 。16kHz，帧长为subplot(211), plot(timit1(14501:14900)%/a/ in dark.subplot(212), plot(timit1(35501:35900)%/s/ in wash.3.3 编写个一用以计算信号

13、短时功率的函数“stpower.m ”，其Matlab代码如下：function Px=stpower (x, N)M = length (x);Px = zeros (M, 1);Px (N)=x (1:N)*x (1:N)/N;For (m = (N+1): M)Px (m) = Px (m-1)+ (x (m) 2-x (m-N) 2)/N;3.4 编写一个用以计算信号短时过零率的函数“stzerocross.m ”，其Matlab代码如下：function Zx = stzerocross(x, N)M=length(x);Zx=zeros(M,1);Zx (N+1)=sum (abs

14、 (sign (x (2:N+1)sign (x (1:N) ) ) )/ (2*N);for (m= (N+2):M)Zx (m) = Zx (m-1) + (abs (sign (x (m) )sign (x (m-1) ) ) abs(sign (x (m-N) ) sign (x (m-N-1 ) ) ) )/ (2*N); end3.5 加载“ digits.mat”语音数据，该数据包含英文数字计算其中的单词“four ”的短时功率和过零率（采样频率为每帧有300 个样点） 。“ 0” 到“ 9”的发音。10kHz ，帧长为30ms ，load digits;N = 300; x=d

15、igits.four1;Px = stpower (x,N);Zx = stzerocross (x,N);plot(Px*1e -5Zxx/2000)3.6 加载“ vowels.mat ”语音数据，该数据包含元音/a/, /i/, /o/, /u/ 的发音。分别画出一帧/a/ 、一帧/i/ 、一帧 /o/和一帧 /u/ 的时域波形（采样频率为10kHz ，帧长为30ms ，每帧有300 个样点） 。load vowelssubplot(221)plot (vowels.a_1(2001:2300)subplot(222)plot (vowels.i_1(2001:2300)subplot(

16、223)plot (vowels.o_1(2001:2300)subplot(224)plot (vowels.u_1(2001:2300)3.7 编写一个利用短时功率和短时过零率这两个参数对语音信号进行浊音分类的函数“ voiunvoi.m ”，其 Matlab 代码如下： function voi = voiunvoi(x,N,Pth,Zth)/清音% Short-time power Px and zero crossing Zx measures. Px = stpower(x,N);Zx = stzerocross(x,N);% Compare estimates with thre

17、shold values.voi = (PxPth*max(Px) & (ZxZth);% Shift the voi-flag N/2 samples to the left (middle sample in window). voi = voi(fix(N/2)+1:length(voi);voi(length(voi)*ones(fix(N/2),1);4、语音信号的频域分析8kHz4.1 加载“ ma1_1 ”语音数据。基于DFT 变换，画出其中一帧数据（采样频率为，帧长为37.5ms ，每帧有300 个样点）的频域波形（对数幅度谱）。load ma1_1;x = ma1_1 (41

18、61:4460);subplot(121), plot(x)N = 1024; k = 0:N/2-1;X = fftshift (fft (x.*hann (length (x),N);subplot(122), plot (k,20*log10 (abs(X(N/2:-1:1),axis (0N/2-1-infinf )4.2 编写一个用以计算信号谐波乘积谱的函数“hpspectrum.m ”，其Matlab代码如下：function HPSx = hpspectrum(x,N,R)k = 1:R:N/2; K = length(k);X = fft (x.*hann(length(x),

19、 N);HPSx = X(k);for (r = R-1:-1:1)HPSx = HPSx.*X (1:r:r*K);end4.3 利用函数“hpspectrum.m ”，来估计一帧浊音的基音周期。代码如下：X = ma1_1(4161:4460);N=1024;R=5;HPSx = hpspectrum (x, N, R);plot (20*log10 (abs (HPSx)4.4 加载“ vowels.mat ”语音数据， 分别画出一帧/i/ 和一帧 /u/（采样频率为10kHz ，帧长为30ms ，每帧有300 个样点）的基于DFT 的对数幅度谱。其Matlab代码如下：load vow

20、elsx = vowels.i_1(2001:2300);N = 1024; k= -N/2:N/2-1;X = fftshift (fft (x.*hann (length(x),N);plot (k,20*log10(abs(X), axis(0 fix (N/2) 0 100)x = vowels.u_1(2001:2300);N= 1024; k = - N/2:N/2-1;X = fftshift (fft (x.*hann(length(x),N);plot (k,20*log10(abs(X), axis(0 fix(N/2) 0 100)4.5 画出一帧清音语音的基于DFT的对

21、数幅度谱。语音数据为15701 15860 个样点（采样频率为8kHz ，帧长为20ms ，每帧有load ma1_1;x = ma1-1 (4161:4460);plot (x)N = 1024; k = - N/2;N/2-1;X = fftshift (fft (x.*hann (length (x),N);plot (k,20*log10 (abs(X), axis (0 fix (N/2)- infinf )ma1_1中的第160 个样点） 。4.6加载“ timit1 ”语音数据。利用函数specgram ，画出该句语音的语谱图。其Matlab代码如下：load timit1;NF

22、FT = 256 ; Fs = 16000; Win = 256; Noverlap = 128;specgram(timit1, NFFT, Fs, Win, Noverlap);4.7 加载“ gliss ”语音数据。画出该句语音的语谱图。load gliss;NFFT = 256 ; Fs = 10000;Win = 256; Noverlap = 128;specgram(gliss.i_2, NFFT, Fs, Win, Noverlap);4.8基于浊音/ 清音分类函数voiunvoi.m和谐波乘积谱函数一个用以计算浊音短时基音周期的函数“stpitch.m ”。其Matlabhp

23、spectrum.m 代码如下：，编写function Fp = stpitch(x,N,Pth,Zth,NFFT,R,Fs)M = length(x);N = 2*fix(N/2);win = hann(N);Fp = zeros(M,1);voi = voiunvoi(x,N,Pth,Zth);invoi = 0;for (m=N:N/2:M)n1 = m-N+1:m-N/2;n2 = m-N/2+1:m;n= n1 n2;if (any(voi(n)if (invoi)Fp(n1) = Fmax;invoi= 0;endelseHPSmax,Fmax = max(abs(hpspectr

24、um(x(n),NFFT,R);if (invoi)Fp(n) = Fp(n) + win*Fmax;elseFp(n) = Fmax*ones(N/2,1) win(N/2+1:N);invoi = 1;endendendwsave = warning; warning(off);Fp = (Fp*Fs/NFFT).*(Fp./Fp);warning(wsave);1. 调用 matlab 中的 wavread 和 wavplay 读入并播放语音文件。2. 调用 fft 命令计算语音的频谱。3. 调用 specgram 命令得到语音的语谱图。4. 调用 plot 命令绘制并观察语音的时频域波

25、形。5. 根据短时功率、短时过零率、谐波乘积谱的计算公式，编写分别用以计算短时 功 率 、 短 时 过 零 率 、 基 音 频 率 的 函 数 ： functionPx=stpower 、 functionZx=stzerocross(x, N)和 function HPSx = hpspectrum(x,N,R)。6. 利用短时功率和短时过零率这两个参数可以对语音信号进行浊音编写实现此分类功能的函数voiunvoi.m。4.9 利用函数stpitch.m ，分别求出两句语音“timit1 ”和“timit2/清音分类。”的基音周期。实验二语音编码一、实验目的熟悉语音基本压缩编码的方法，观察语

26、音压缩效果，加深对语音线性预测编码（ LPC ）的理解。二、实验内容1、编写并调试语音LPC 参数提取程序。2、编写并调试语音基音周期提取程序。3、编写并调试语音LPC 合成程序。三、实验原理语音信号中含有大量的冗余信息，采用各种信源编码技术减除语音信号的冗余度，并充分利用人耳的听觉掩蔽效应，就可以将其编码速率压缩很多倍，而仍能提供可懂语音。LPC 声码器是一种比较简单实用的语音压缩方法，其基本原理是：根据语音生成模型，将语音看作激励源通过一个线性时不变系统产生的输出，利用线性预测分析对声道参数进行估值，将求得的线性预测系数，结合基音周期等少量参数进行传输，就可以在接收端利用合成滤波器重构语音

27、信号。线性预测系数的估计方法为：假设语音的当前样值可以用过去的M 个语音样值来进行预测MaM x n Mai x n ix n a1 x n 1 a2 x n 2i 1式中ai 即为线性预测系数。实际值和预测值之间的均方误差可表示为M22 nEx nai x n inni 1要求均方误差总和最小，将E 关于 ai 的偏导数设置为零，可以得到Mx nkx nai x ni0ni1通过采用自相关法、协方差法或格形法求解该方程，即可得到最优的ai。四、实验方法及程序1. 调用 xcorr 命令计算一帧语音的自相关函数。2. 调用 toeplitz命令形成该帧语音的自相关矩阵。3. 调用 durbin

28、 命令，采用杜宾递推算法计算该帧语音的线性预测系数。4. 编写 lpcauto.m 函数，求取一句语音信号的线性预测系数及预测残差。选择设当的窗函数对语音信号进行分幀。5. 编写 lpcpitch 函数，由残差信号计算该句语音的基音周期。6. 编写 lpcgain 函数，由预测残差能量，求出该句语音的增益。7. 编写 lpcsyn 函数，由该句语音的基音周期、预测残差能量和增益进行LPC 合成。五、实验结果与分析1. 如 何对 全极点 模型的线 性预测 参数进 行提取？ 加窗等 预处理 对提取结 果有何影响？自相关法、协方差法求解线性预测方程组、格型法。协方差法是先把计算均方误差的间隔 N固定

29、下来，再研究它对计算 n（ i,1 ）的其中 N 是基音周期值。添加汉明窗和汉宁窗可以得到效果不错的波形，很明显两个波形相似，矩形窗的波形相对不是很平稳。通过改变帧长、帧移后波形更加的清晰，特征更加的明显，更方便观察。2. 如何实现线性预测参数的各种表现方式之间的转换？对于语音信号，确定了各线性预测系数后，根据H （ z）可得其频率响应的估值即 LPC 谱。3. 使用预测残差求基音周期有何优点？该方法能进一步去除声道共振峰的影响，简化清浊音判决过程，提高判决精度，有效降低发生基音倍频和半频的错误。4. 如何根据线性预测系数求得LPC 频谱？对于语音信号，确定了各线性预测系数后，根据H （ z）

30、可得其频率响应的估值即 LPC 谱。5. LPC 频谱与基于 DFT 求得的语音对数幅度谱相比有何不同？LPC 谱优点在于可以很好地表示共振峰结构而不出现额外的蜂起和起伏，基本不含有噪声，其频率灵敏度与人耳不相匹配； DET 得到的对数幅度谱受基频谐波的影响，最大值只出现在谐波频率上，共振峰测定误差较大。6. LPC 合成语音的质量如何？有何改进措施？LPC 合成技术的优点是简单直观。其合成过程实质上只是一种简单的解码和拼接过程。另外，由于波形拼接技术的合成基元是语音的波形数据，保存了语音的全部信息，因而对于单个合成基元来说能够获得很高的自然度。但是，由于自然语流中的语音和孤立状况下的语音有着

31、极大的区别，如果只是简单地把各个孤立的语音生硬地拼接在一起，其整个语流的质量势必是不太理想的。而LPC 技术从本质上来说只是一种录音+ 重放，对于合成整个连续语流LPC 合成技术的效果是不理想的。因此， LPC 合成技术必须和其他技术相结合，才能明显改善LPC 合成的质量。实验过程记录与结果分析1、计算短时自相关函数1.1分别计算正弦信号和白噪声的短时自相关函数，估计正弦信号的基音周期。其 Matlab 代码如下：x=sin(2*pi*0.01*(0:499);r,eta=xcorr(x,100,unbiased);stem(eta,r);w=randn(500,1);r,eta=xcorr(

32、w,100,unbiased);stem(eta,r);1.2计算正弦信号叠加白噪声的短时自相关函数，试估计正弦信号的基音周期。其 Matlab 代码如下：x=sin(2*pi*0.01*(0:499);w=randn(500,1);x1=x+w;r,eta=xcorr(x1,100,unbiased);stem(eta,r);长为1.3分别画出一帧浊音和一帧清音的语音自相关波形（采样频率为25.6ms ，每帧有256 个样点） ，估计浊音的基音周期。load digits; x=digits.three1;m=2756; N=256; n=m-N+1:m;r,eta=xcorr(x(n),2

33、50,unbiased);plot(eta,r);10kHz，帧m=500; N=256; n=m-N+1:m;r,eta=xcorr(x(n),250,unbiased);plot(eta,r);2、求取 LPC系数2.1加载“ digits.three1”语音数据。基于自相关法，求出其中一帧数据（采样频率为10kHz ，帧长为25.6ms ，每帧有256 个样点）的14 阶 LPC 系数。画出相应的 LPC 谱。load digits; x=digits.three1;m=2756; N=256; n=m-N+1:m; M=14;r,eta=xcorr(x(n),250,unbiased)

34、;Rx=toeplitz(r(M+1:2*M);rx=r(M+2:2*M+1);a=Rxrx;NFFT=1024; k=1:NFFT/2;X=fft(x(n).*hann(N),NFFT);Theta=1./fft(1;-a,NFFT);plot(k,20*log10(abs(353*Theta(k) X(k)axis(0 NFFT/2inf inf)2.2 编写一个用以实现杜宾递推算法的函数“durbin.m ”，其 Matlab代码如下：function a,xi,kappa=durbin(r,M)kappa = zeros (M,1);a = zeros (M, 1);xi=r(1);

35、zeros(M,1);for (j=1:M)kappa(j)=r(j+1)-a(1:j-1)*r(j:-1:2)/xi(j);a(j)=kappa(j);a(1:j-1)=a(1:j-1)-kappa(j)*a(j-1:-1:1);xi(j+1)=a(1:j-1)-kappa(j)*a(j-1:-1:1);xi(j+1)=xi(j)*(1-kappa(j)2);end2.3 利用函数“durbin.m ”，计算2.1 中语音数据帧的14 阶 LPC 系数，并与2.1中的结果进行比较：x=digits.three1;m=2756; N=256; n=m-N+1:m; M=14;r,eta=xco

36、rr(x(n),250,unbiased);aLD,xi,kappa=durbin(r(M+1:2*M+1),M);a, aLD, norm(a-aLD)a =0.35120.20720.38080.31630.0522-0.1639-0.2903-0.2898-0.12220.19600.34920.30240.1713-0.8258aLD =0.35120.20720.38080.31630.0522-0.1639-0.2903-0.2898-0.12220.19600.34920.30240.1713-0.8258ans =1.0047e-142.4 编写一个用以实现反射系数转换为LPC

37、 系数的函数“rf2lpc.m ”，其Matlab代码如下：function a = rf2lpc(kappa)M = length(kappa);a = zeros(M,1);for (j=1:M)a(j) = kappa(j);a(1:j-1) = a(1:j-1) - kappa(j)*a(j-1:-1:1);end2.5编写一个用以实现LPC 系数转换为反射系数的函数“lpc2rf.m”，其 Matlab代码如下：function kappa = lpc2rf(a)M = length(a);kappa = zeros(M,1);for (j=M:-1:1)kappa(j) = a(j

38、);a(1:j-1) = (a(1:j-1) + a(j)*a(j-1:-1:1)/(1 - kappa(j)2);end2.6 使用函数“rf2lpc.m ”和“ lpc2rf.m ”，检验反射系数和LPC 系数相互转换结果。norm(kappa-lpc2rf(aLD)norm(aLD-rf2lpc(kappa)ans =1.4760e-14ans =03、语音信号的逆滤波3.1 利用 2.3 中求出的一帧语音数据的 LPC 系数，构造逆滤波器，并画出该帧语音信号的残差波形。ehat=filter(1;-1,1,x(n);plot(x(n) ehat)3.2 编写一个用以实现语音信号LPC

39、分析的函数“lpcauto ”，其Matlab代码如下：function ar,xi,e,m = lpcauto(x,M,win,Olap)Nx = length(x);N = length(win);if (N = 1)N = win;win = ones(N,1);endF = fix(Nx-Olap)/(N-Olap);ar = zeros(M+1,F);xi = zeros(M+1,F); e = zeros(Nx,1); m = zeros(F,1);n = 1:N;n1 = 1:Olap;n2 = N-Olap+1:N;n3 = Olap+1:N;win1 = win(n1)./(win(n1)+win(n2)+