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

-.z.邮电大学通达学院实验报告实验名称语音时频参数的提取和应用实验语音编码实验课程名称信息处理技术专业综合实验班级**姓名开课时间2016/2017学年,第二学期实验一语音时频参数的提取和应用一、实验目的熟悉根本的语音时频参数提取方法，观察并比拟各类参数采用不同提取方法的效果，了解这些时频参数在语音处理中的应用。二、实验容1、编写并调试语音时域参数提取程序，包括短时功率、短时过零率等函数。2、编写并调试语音频域参数提取程序，包括频谱、语谱图、基音频率、共振峰等。3、利用提取出的时频参数对语音信号进展静/清/浊音的划分。三、实验原理语音信号具有短时平稳性，需对其进展分段〔帧〕处理〔帧长一般取10至30ms之间〕。根据鼓励方式的不同，语音信号可以分成浊音和清音两大类：浊音：鼓励源可模拟成准周期性的脉冲串。该周期称为基音周期，其倒数就是基音频率，基音频率一般在几百赫兹以。清音：鼓励源可模拟成随机白噪声。简化的语音生成数学模型如下图：对语音时\频波形进展绘制与观察，有助于了解语音的清/浊音等特性。对基音周期等声门参数和共振峰等声道参数进展提取与分析，是语音处理与应用的重要容。对于时变语音信号，可以计算它的短时功率和短时过零率。第m帧〔帧长为N个样点〕的功率计算公式为：过零率计算公式为：对于浊音语音，可以利用其频谱具有丰富的谐波分量的特点，求出其谐波乘积谱，计算公式为：式中，R一般取为5。在谐波乘积谱中，基频分量变得很大，更易于估计基音周期。四、实验方法及程序1.调用matlab中的wavread和wavplay读入并播放语音文件。2.调用fft命令计算语音的频谱。3.调用specgram命令得到语音的语谱图。4.调用plot命令绘制并观察语音的时\频域波形。5.根据短时功率、短时过零率、谐波乘积谱的计算公式，编写分别用以计算短时功率、短时过零率、基音频率的函数：functionP*=stpower、functionZ*=stzerocross(*,N)和functionHPS*=hpspectrum(*,N,R)。6.利用短时功率和短时过零率这两个参数可以对语音信号进展浊音/清音分类。编写实现此分类功能的函数voiunvoi.m。五、实验结果与分析1.清音和浊音的短时功率、短时过零率各有何特点？清音段的能量一般比浊音段的小得多。浊音时能量集中于较低频率段，具有较低的过零率，而清音时能量集中于较高频率段，具有较高的过零率。2.如何利用短时功率和短时过零率这两个参数对语音信号进展浊音/清音的分类？发浊音时，由于声门波引起了谱的高频跌落。所以其语音能量集中在3kHz以下；发清音时，多数能量出现在高频。高频率意味着高的平均过零率，低频率意味着低的平均过零率，一般可认为浊音具有较低的平均过零率，清音具有较高的平均过零率。3.如何根据时域波形估计各帧元音的基音周期？看图可知四帧语音的周期均为0.01s，第一帧0~100Hz有10个峰值点，基频约为10Hz，第二帧0~100Hz有5个峰值点，基频约为20Hz，第三帧0~100Hz有5个峰值点，基频约为20Hz，第四帧0~100Hz有4个峰值点，基频约为25Hz。4.如何根据基于DFT的对数幅度谱估计*幀浊音的基音周期？0~1172为一个基因周期。5.如何根据基于DFT的对数幅度谱，估计出共振峰频率？有7个峰点，因此基频约为167。6.时域对语音信号进展加窗，反映在频域，其窗谱对基于DFT的对数幅度谱有何影响？如何估计出窗谱的主瓣宽度？因为时间窗幅度的傅立叶变换为脉冲状态,则频谱偏差很小,当定常过程*(n)的平均值β为零且N很大时E[TN(ω)]=f(ω)(12)所以加窗的周期图TN(ω)是频谱密度函数的无偏估计。固有TN(ω)=[1+ξ(ω)]f(ω)。如果原始信号的频谱成份与FFT中的谱线完全一致，这种情况下采样数据的长度为信号周期的整数倍，频谱中只有主瓣。7.如何根据窄带、宽带语谱图，提取出语音的基音和共振峰轨迹？在基音提取中，广泛采用语音波形或误差信号波形的低通滤波。提取共振峰特性最简便的手段是使用语谱仪。提取共振峰还有倒谱法、LPC分析法等更为有效、准确的方法。〔做完试验后，记住要答复这些问题〕实验过程记录与结果分析1、进入matlab1.1在Windows环境中，创立一个名为speech的文件夹，存放语音数据和与实验相关的Matlab文件。如：E:\speech。用于实验的语音数据〔.mat〕及相应说明〔.t*t〕包括：digits.mat %英文数字“0”到“9digits.t*tgliss.mat%两个包含/i/的滑音gliss.t*tletters.mat%英文字母表中26个字母的发音letters.t*tma1.t*tma1-1.mat%一句连续语音ma1-2.mat%一句连续语音ma1-3.mat%一句连续语音ma1-4.mat%一句连续语音timit.t*ttimit1.mat%一句连续语音timit2.mat%一句连续语音timit3.mat%一句连续语音timit4.mat%一句连续语音vowels.mat%元音/a/,/i/,/o/,/u/的发音vowels.t*twords.mat%十个孤立字的发音words.t*t1.2在Windows桌面上，双击Matlab的图标，进入Matlab的工作环境。在命令窗口>>提示符后，键入cdE:\speech指令，进入实验文件夹路径。1.3在命令窗口键入clearall指令，去除工作空间的所有变量。2、加载语音数据，熟悉matlab命令2.1加载“timit1.mat〞语音数据，命令为：loadtimit12.2用命令“who〞和“whos〞列出当前工作空间中的变量。变量“timit1”Yourvariablesare:timit1NameSizeBytesClassAttributestimit155911*1447288double2.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-0.02000.20210.40120.1869-0.1756-0.19020.08230.23630.17120.0637-0.0209-0.0683-0.03940.00810.01050.03110.06962.5Matlab中，符号“’〞可以实现转置运算，运行以下命令进展观察：size(timit1’〔粘贴实验结果〕ans=1559112.6使用函数ma*和min可以分别得到信号的最大值和最小值：ma*(timit1)min(timit1)〔粘贴实验结果〕ans=1ans=-0.80683、语音信号的时域分析3.1使用以下命令，画出语音信号timit1〔采样频率为16kHz，长度为3.5秒〕的时域波形。plot(timit1)title(‘timit1’*label(‘TimeInde*,{\itn}’)ylabel(‘Amplitude’)a*is([0length(timit1)-11])grid3.2分别画出一帧浊音和一帧清音的语音时域波形〔采样频率为16kHz，帧长为25ms，每帧有400个样点〕。subplot(211),plot(timit1(14501:14900))%/a/indark.subplot(212),plot(timit1(35501:35900))%/s/inwash.3.3编写个一用以计算信号短时功率的函数“stpower.m〞，其Matlab代码如下：functionP*=stpower(*,N)M=length(*);P*=zeros(M,1);P*(N)=*(1:N)’**(1:N)/N;For(m=(N+1):M)P*(m)=P*(m-1)+(*(m)^2-*(m-N)^2)/N;3.4编写一个用以计算信号短时过零率的函数“stzerocross.m〞，其Matlab代码如下：functionZ*=stzerocross(*,N)M=length(*);Z*=zeros(M,1);Z*(N+1)=sum(abs(sign(*(2:N+1))–sign(*(1:N))))/(2*N);for(m=(N+2):M)Z*(m)=Z*(m-1)+(abs(sign(*(m))–sign(*(m-1)))…－abs(sign(*(m-N))–sign(*(m-N-1))))/(2*N);end3.5加载“digits.mat〞语音数据，该数据包含英文数字“0”到“9”的发音。计算其中的单词loaddigits;N=300;*=digits.four1;P*=stpower(*,N);Z*=stzerocross(*,N);plot([P**1e-5Z**/2000])3.6加载“vowels.mat〞语音数据，该数据包含元音/a/,/i/,/o/,/u/的发音。分别画出一帧/a/、一帧/i/、一帧/o/和一帧/u/的时域波形〔采样频率为10kHz，帧长为30ms，每帧有300个样点〕。loadvowelssubplot(221)plot(vowels.a_1(2001:2300))subplot(222)plot(vowels.i_1(2001:2300))subplot(223)plot(vowels.o_1(2001:2300))subplot(224)plot(vowels.u_1(2001:2300))3.7编写一个利用短时功率和短时过零率这两个参数对语音信号进展浊音/清音分类的函数“voiunvoi.m〞，其Matlab代码如下：functionvoi=voiunvoi(*,N,Pth,Zth)%Short-timepowerP*andzerocrossingZ*measures.P*=stpower(*,N);Z*=stzerocross(*,N);%pareestimateswiththresholdvalues.voi=(P*>Pth*ma*(P*))&(Z*<Zth);%Shiftthevoi-flagN/2samplestotheleft(middlesampleinwindow).voi=[voi(fi*(N/2)+1:length(voi));voi(length(voi))*ones(fi*(N/2),1)];4、语音信号的频域分析4.1加载“ma1_1”loadma1_1;*=ma1_1(4161:4460);subplot(121),plot(*)N=1024;k=0:N/2-1;*=fftshift(fft(*.*hann(length(*)),N));subplot(122),plot(k,20*log10(abs(*(N/2:-1:1)))),a*is([0N/2-1-infinf])4.2编写一个用以计算信号谐波乘积谱的函数“hpspectrum.m〞，其Matlab代码如下：functionHPS*=hpspectrum(*,N,R)k=1:R:N/2;K=length(k);*=fft(*.*hann(length(*)),N);HPS*=*(k);for(r=R-1:-1:1)HPS*=HPS*.**(1:r:r*K);end4.3利用函数“hpspectrum.m〞，来估计一帧浊音的基音周期。代码如下：*=ma1_1(4161:4460);N=1024;R=5;HPS*=hpspectrum(*,N,R);plot(20*log10(abs(HPS*)))4.4加载“vowels.mat〞语音数据，分别画出一帧/i/和一帧/u/〔采样频率为10kHz，帧长为30ms，每帧有300个样点〕的基于DFT的对数幅度谱。其Matlab代码如下：loadvowels*=vowels.i_1(2001:2300);N=1024;k=-N/2:N/2-1;*=fftshift(fft(*.*hann(length(*)),N));plot(k,20*log10(abs(*))),a*is([0fi*(N/2)0100])*=vowels.u_1(2001:2300);N=1024;k=-N/2:N/2-1;*=fftshift(fft(*.*hann(length(*)),N));plot(k,20*log10(abs(*))),a*is([0fi*(N/2)0100])4.5画出一帧清音语音的基于DFT的对数幅度谱。语音数据为ma1_1中的第15701－15860个样点〔采样频率为8kHz，帧长为20ms，每帧有160个样点〕。loadma1_1;*=ma1-1(4161:4460);plot(*)N=1024;k=-N/2;N/2-1;*=fftshift(fft(*.*hann(length(*)),N));plot(k,20*log10(abs(*))),a*is([0fi*(N/2)-infinf])4.6加载“timit1”loadtimit1;NFFT=256;Fs=16000;Win=256;Noverlap=128;specgram(timit1,NFFT,Fs,Win,Noverlap);4.7加载“gliss〞语音数据。画出该句语音的语谱图。loadgliss;NFFT=256;Fs=10000;Win=256;Noverlap=128;specgram(gliss.i_2,NFFT,Fs,Win,Noverlap);4.8基于浊音/清音分类函数voiunvoi.m和谐波乘积谱函数hpspectrum.m，编写一个用以计算浊音短时基音周期的函数“stpitch.m〞。其Matlab代码如下：functionFp=stpitch(*,N,Pth,Zth,NFFT,R,Fs)M=length(*);N=2*fi*(N/2);win=hann(N);Fp=zeros(M,1);voi=voiunvoi(*,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=[n1n2];if(any(~voi(n)))if(invoi)Fp(n1)=Fma*;invoi=0;endelse[HPSma*,Fma*]=ma*(abs(hpspectrum(*(n),NFFT,R)));if(invoi)Fp(n)=Fp(n)+win*Fma*;elseFp(n)=Fma**[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命令绘制并观察语音的时\频域波形。5.根据短时功率、短时过零率、谐波乘积谱的计算公式，编写分别用以计算短时功率、短时过零率、基音频率的函数：functionP*=stpower、functionZ*=stzerocross(*,N)和functionHPS*=hpspectrum(*,N,R)。6.利用短时功率和短时过零率这两个参数可以对语音信号进展浊音/清音分类。编写实现此分类功能的函数voiunvoi.m。4.9利用函数stpitch.m，分别求出两句语音“timit1”和“timit2实验二语音编码一、实验目的熟悉语音根本压缩编码的方法，观察语音压缩效果，加深对语音线性预测编码〔LPC〕的理解。二、实验容1、编写并调试语音LPC参数提取程序。2、编写并调试语音基音周期提取程序。3、编写并调试语音LPC合成程序。三、实验原理语音信号中含有大量的冗余信息，采用各种信源编码技术减除语音信号的冗余度，并充分利用人耳的听觉掩蔽效应，就可以将其编码速率压缩很多倍，而仍能提供可懂语音。LPC声码器是一种比拟简单实用的语音压缩方法，其根本原理是：根据语音生成模型，将语音看作鼓励源通过一个线性时不变系统产生的输出，利用线性预测分析对声道参数进展估值，将求得的线性预测系数，结合基音周期等少量参数进展传输，就可以在接收端利用合成滤波器重构语音信号。线性预测系数的估计方法为：假设语音的当前样值可以用过去的个语音样值来进展预测式中即为线性预测系数。实际值和预测值之间的均方误差可表示为要求均方误差总和最小，将关于的偏导数设置为零，可以得到通过采用自相关法、协方差法或格形法求解该方程，即可得到最优的。四、实验方法及程序1.调用*corr命令计算一帧语音的自相关函数。2.调用toeplitz命令形成该帧语音的自相关矩阵。3.调用durbin命令，采用杜宾递推算法计算该帧语音的线性预测系数。4.编写lpcauto.m函数，求取一句语音信号的线性预测系数及预测残差。选择设当的窗函数对语音信号进展分幀。5.编写lpcpitch函数，由残差信号计算该句语音的基音周期。6.编写lpcgain函数，由预测残差能量，求出该句语音的增益。7.编写lpcsyn函数，由该句语音的基音周期、预测残差能量和增益进展LPC合成。五、实验结果与分析1.如何对全极点模型的线性预测参数进展提取？加窗等预处理对提取结果有何影响？自相关法、协方差法求解线性预测方程组、格型法。协方差法是先把计算均方误差的间隔N固定下来，再研究它对计算n〔i,1〕的其中N是基音周期值。添加汉明窗和汉宁窗可以得到效果不错的波形，很明显两个波形相似，矩形窗的波形相对不是很平稳。通过改变帧长、帧移后波形更加的清晰，特征更加的明显，更方便观察。2.如何实现线性预测参数的各种表现方式之间的转换？对于语音信号，确定了各线性预测系数后，根据H〔z〕可得其频率响应的估值即LPC谱。3.使用预测残差求基音周期有何优点？该方法能进一步去除声道共振峰的影响，简化清浊音判决过程，提高判决精度，有效降低发生基音倍频和半频的错误。4.如何根据线性预测系数求得LPC频谱？对于语音信号，确定了各线性预测系数后，根据H〔z〕可得其频率响应的估值即LPC谱。5.LPC频谱与基于DFT求得的语音对数幅度谱相比有何不同？LPC谱优点在于可以很好地表示共振峰构造而不出现额外的蜂起和起伏，根本不含有噪声，其频率灵敏度与人耳不相匹配；DET得到的对数幅度谱受基频谐波的影响，最大值只出现在谐波频率上，共振峰测定误差较大。6.LPC合成语音的质量如何？有何改良措施？LPC合成技术的优点是简单直观。其合成过程实质上只是一种简单的解码和拼接过程。另外，由于波形拼接技术的合成基元是语音的波形数据，保存了语音的全部信息，因而对于单个合成基元来说能够获得很高的自然度。但是，由于自然语流中的语音和孤立状况下的语音有着极大的区别，如果只是简单地把各个孤立的语音生硬地拼接在一起，其整个语流的质量势必是不太理想的。而LPC技术从本质上来说只是一种录音+重放，对于合成整个连续语流LPC合成技术的效果是不理想的。因此，LPC合成技术必须和其他技术相结合，才能明显改善LPC合成的质量。实验过程记录与结果分析1、计算短时自相关函数1.1分别计算正弦信号和白噪声的短时自相关函数，估计正弦信号的基音周期。其Matlab代码如下：*=sin(2*pi*0.01*(0:499)');[r,eta]=*corr(*,100,'unbiased');stem(eta,r);w=randn(500,1);[r,eta]=*corr(w,100,'unbiased');stem(eta,r);1.2计算正弦信号叠加白噪声的短时自相关函数，试估计正弦信号的基音周期。其Matlab代码如下：*=sin(2*pi*0.01*(0:499)');w=randn(500,1);*1=*+w;[r,eta]=*corr(*1,100,'unbiased');stem(eta,r);1.3分别画出一帧浊音和一帧清音的语音自相关波形〔采样频率为10kHz，帧长为25.6ms，每帧有256个样点〕，估计浊音的基音周期。loaddigits;*=digits.three1;m=2756;N=256;n=m-N+1:m;[r,eta]=*corr(*(n),250,'unbiased');plot(eta,r);m=500;N=256;n=m-N+1:m;[r,eta]=*corr(*(n),250,'unbiased');plot(eta,r);2、求取LPC系数2.1加载“digits.three1”loaddigits;*=digits.three1;m=2756;N=256;n=m-N+1:m;M=14;[r,eta]=*corr(*(n),250,'unbiased');R*=toeplitz(r(M+1:2*M));r*=r(M+2:2*M+1);a=R*\r*;NFFT=1024;k=1:NFFT/2;*=fft(*(n).*hann(N),NFFT);Theta=1./fft([1;-a],NFFT);plot(k,20*log10(abs([353*Theta(k)*(k)])))a*is([0NFFT/2–infinf])2.2编写一个用以实现杜宾递推算法的函数“durbin.m〞，其Matlab代码如下：function[a,*i,kappa]=durbin(r,M)kappa=zeros(M,1);a=zeros(M,1);*i=[r(1);zeros(M,1)];for(j=1:M)kappa(j)=r(j+1)-a(1:j-1)'*r(j:-1:2)/*i(j);a(j)=kappa(j);a(1:j-1)=a(1:j-1)-kappa(j)*a(j-1:-1:1);*i(j+1)=a(1:j-1)-kappa(j)*a(j-1:-1:1);*i(j+1)=*i(j)*(1-kappa(j)^2);end2.3利用函数“durbin.m〞，计算2.1中语音数据帧的14阶LPC系数，并与2.1中的结果进展比拟：*=digits.three1;m=2756;N=256;n=m-N+1:m;M=14;[r,eta]=*corr(*(n),250,'unbiased');[aLD,*i,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系数的函数“rf2lpc.m〞，其Matlab代码如下：functiona=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代码如下：functionkappa=lpc2rf(a)M=length(a);kappa=zeros(M,1);for(j=M:-1:1)kappa(j)=a(j);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,*(n));plot([*(n)ehat])3.2编写一个用以实现语音信号LPC分析的函数“lpcauto〞，其Matlab代码如下：function[ar,*i,e,m]=lpcauto(*,M,win,Olap)N*=length(*);N=length(win);if(N==1)N=win;win=ones(N,1);endF=fi*((N*-Olap)/(N-Olap));ar=zeros(M+1,F);*i=zeros(M+1,F);e=zeros(N*,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)+eps);win2=win(n2)./(win(n1)+win(n2)+eps);for(f=1:F)[r,eta]=*corr(*(n).*win,M,'biased');[a,*i(:,f)