语音信号处理实验

实验一语音信号的预处理1实验目的通过Matlab编程掌握语音信号的预处理方法，包括对信号进行分帧、预加重，加窗处理。2实验原理由于语音信号从整体上来看是一个非平稳过程，但是在一个短的时间内，其特性保持相对不变，所以语音信号具有短时平稳性，对语音信号的分析和处理必须建立在“短时”的基础上，将信号分为一段一段来分析其特征参数。人发音时存在口唇的辐射效应，口唇的辐射模型相当于一阶高通滤波器，所以在对实际信号进行分析处理时，常用“预加重技术”，目的提升信号的高频部分，使信号的频谱更加平坦，方便信号的分析。3实验过程1)读语音数wavread2)听语音sound3)写语音wavwrite4）对语音信号进行分帧处理5）对语音信号进行预加重1()10.96Hzz−=−6）加汉明窗0.420.5*cos(2/)0.08*cos(2/)0()0kMkMkMwkππ−+≤≤⎧=⎨⎩其他MkMkkw≤≤+−=1)],12cos(1[5.0)(π4实验结果实验一程序所用函数wavread、sound、wavwrite、length、fix、y=x(m:n)、subplot、plot、figure、xlabel、ylabel、title选取c:\windows\media目录下.wav文件作为读取的声音文件。读[x,fs,bits]=wavread('c:\windows\media\dig.wav');听sound(x,fs)写wavwrite(x,fs,bits,'c:\windows\media\dig.wav')分帧处理流程1读语音数据。2求语音长度。3确定帧长和帧移。4确定可以分多少帧5for循环实现各帧信号的取出。先确定各帧信号的起点和终点坐标，然后利用矩阵操作函数将各帧信号取出。[x,fs,bits]=wavread('c:\windows\media\dig.wav');x=x(:,1);x=x';len=length(x);N=256;M=128;Fn=fix((len-N)/M+1);y=[];fori=1:Fndown=1+(i-1)*M;up=down+N-1;temp=x(down:up);y=[y;temp];end加窗取出一帧信号和窗函数相乘预加重Z域1()10.96Hzz−=−y=filter([1-0.96],1,x);时域表达：y(k)=x(k)-x(k-1)程序：[x,fs,bits]=wavread('c:\windows\media\dig.wav');x=x(:,1);y=x(10101:10612);figure(1)plot(y)L=length(y);k=1:L;w=0.42-0.5*cos(2*pi*k/L)+0.08*cos(4*pi*k/L);z=y.*w’;figure(2)plot(z)forr=2:Lp(r)=y(r)-y(r-1);endfigure(3)plot(p)实验二语音信号的时域分析1实验目的通过Matlab编程实现语音信号的时域波形图，并观察清音、浊音信号的时域特点。掌握语音信号的时域分析技术，如短时平均能量、短时平均幅度、短时平均过零率分析、短时平均自相关、短时平均幅度差。2实验原理语音信号的时域分析就是分析和提取语音信号的时域参数。时域分析通常用于最基本的参数分析及应用，此方法表示语音信号比较直观，物理意义明确，实现起来比较简单，运算量少。3实验过程1)观察浊音信号波形图2)观察清音信号波形图3)计算语音信号的短时能量、短时平均幅度并画图120()NnnmExm−==∑10|()|NnnmMxm−==∑4)计算信号的短时平均过零率并画图101{|sgn[()]sgn[(1)]|}2NnnmZnxmxm−==−−∑5)计算语音信号的短时自相关函数并画图10()()()NknnnmRkxmxmk−−==+∑6)计算语音信号的短时平均幅度差函数并画图10()|()()|NknnnmFkxmxmk−−==−+∑4实验结果1实验二程序%短时能量[x,fs,bits]=wavread('c:\WINDOWS\Media\chimes.wav');x=x(:,1);x=x';len=length(x);N=256;M=128;Fn=fix((len-N)/M+1);y=[];fori=1:Fndown=1+(i-1)*M;up=down+N-1;temp=x(down:up);y=[y;temp];end%K=100;E=[];fori=1:Fntemp=sum(y(i,:).^2,2);E=[E,temp];endplot(E)%短时平均幅度差[x,fs,bits]=wavread('c:\WINDOWS\Media\chimes.wav');x=x(:,1);x=x';len=length(x);N=256;M=128;Fn=fix((len-N)/M+1);y=[];fori=1:Fndown=1+(i-1)*M;up=down+N-1;temp=x(down:up);y=[y;temp];end%K=100;M=[];fori=1:Fntemp=sum(abs(y(i,:)),2)2M=[M,temp];endplot(M)%短时过零率[x,fs,bits]=wavread('c:\WINDOWS\Media\chimes.wav');x=x(:,1);x=x';len=length(x);N=256;M=128;Fn=fix((len-N)/M+1);y=[];fori=1:Fndown=1+(i-1)*M;up=down+N-1;temp=x(down:up);y=[y;temp];end%K=100;z=[];fori=1:Fntemp=sum(abs(sign(y(i,2:N))-sign(y(i,1:N-1))),2);z=[z,temp];endplot(z)%短时自相关函数[x,fs,bits]=wavread('c:\WINDOWS\Media\chimes.wav');x=x(:,1);x=x';len=length(x);N=256;M=128;Fn=fix((len-N)/M+1);y=[];fori=1:Fndown=1+(i-1)*M;up=down+N-1;temp=x(down:up);y=[y;temp];end3K=100;fori=1:Fnforj=0:KR(i,j+1)=sum(y(i,1:N-j).*y(i,j+1:N),2);endendplot(R(111,:));%短时平均幅度差[x,fs,bits]=wavread('c:\WINDOWS\Media\chimes.wav');x=x(:,1);x=x';len=length(x);N=256;M=128;Fn=fix((len-N)/M+1);y=[];fori=1:Fndown=1+(i-1)*M;up=down+N-1;temp=x(down:up);y=[y;temp];end%K=100;fori=1:Fnforj=0:KF(i,j+1)=sum(abs(y(i,1:N-j)-y(i,j+1:N)),2);endendplot(F(111,:));实验三语音信号的mel频率倒谱参数1实验目的通过Matlab编程掌握语音信号的mel频率倒谱参数的求解方法。2实验原理人耳听到声音的高低与声音的频率成对数关系，即：Mel(f)=2595lg(1+f/700)，实际频率f的单位是Hz.根据Zwicker的工作，临界频率带宽随着频率的变化而变化，并与Mel频率的增长一致。类似与临界带的划分，可以将语音频率划分成一系列三角形的滤波器序列，如下图所示.取每个三角形的滤波器频率带宽内所有信号幅度加权和作为某个带通滤波器的输出，然后对所有滤波器输出做对数运算，再进一步做离散余弦变换得到MFCC。具体步骤如下：（1）三角滤波器的输出则为此频率带宽内所有信号幅度谱加权和。()()()()()()()|()||()|()()()()clhlnnkolkclkolhlkYlXkXkclolhlcl==−−=+−−∑∑l=1,2,....,40（2）对所有滤波器输出作对数运算ln(())Yll=1,2,....,40（3）作离散余弦变换（DCT）得到Mel频率倒谱参数(MFCC)。2411ln(())cos[()]224ilCYlilπ==−∑i=1,2,...,P，P为MFCC参数的阶数，取P＝16。3实验过程4实验结果实验三程序[x,fs,bits]=wavread('c:\WINDOWS\Media\chimes.wav');x=x(:,1);x=x';len=length(x);N=256;M=128;Fn=fix((len-N)/M+1);y=[];fori=1:Fndown=1+(i-1)*M;up=down+N-1;temp=x(down:up);temp=temp.*hamming(N)';%每帧加hamming窗y=[y;temp];endL=24;%滤波器个数R=16;%mel参数的个数k=0:N/2;f=fs/N*k;%每点的频率值mel=2595*log(1+f/700);%mel频率的变换%melm=max(mel)melm=2595*log(1+fs/1400);%mel频率的最高频率r=0:L+1;%整个频率分为r个临界点tri=melm/(L+1)*r;%每个mel临界频率的频率值s=[];forj=1:Fntemp1=y(j,:);%提取一帧信号p=abs(fft(temp1));%每帧信号的幅度谱forl=1:Ltri1=[tri(l),tri(l+1),tri(l+2)];%三角滤波器的下截止频率，中心频率，上截止频率low=find((mel=tri1(1))&(mel=tri1(2)));%三角滤波器下截止频率high=find((mel=tri1(2))&(mel=tri1(3)));%三角滤波器上截止频率w=[(mel(low)-tri1(1))/(tri1(2)-tri1(1)),(tri1(3)-mel(high))/(tri1(3)-tri1(2))];%三角形滤波器m(l)=sum(w.*p([low,high]),2);%给每一段临界频率加三角窗endl=1:L;forq=1:Rc(q)=sqrt(2/N)*sum(log(m).*cos((l-0.5)*q*pi/L),2);%离散余弦变换得到MFCCends=[s;c];endplot(s(6,:));%输出某一帧的MFCC参数实验四语音信号的LPC倒谱参数的提取1实验目的通过Matlab编程掌握语音信号的LPC倒谱参数的求解方法。2实验原理线性预测分析的基本思想是：用过去P个样点的值来预测现在或未来的样点值：预测误差为：理论上通常采用预测系数{ak}的取值使均方误差E[e2(n)](是{ak}的函数)最小的准则。]))()([()]([122∑=−−=pkkknxanxEneE0)]([2=∂∂kaneEpk≤≤10)]()([0)]()([2)]([2=−=−−=∂∂jnxneEknxneEaneEk预测误差与过去的样点不相关[]0)()()()()]()([)(])()([0)()(111=−−=⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡−⋅−−−=⎭⎬⎫⎩⎨⎧−⋅−−⇒=−∑∑∑===pkkpkkpkkkjRajRjnxknxajnxnxEjnxknxanxEjnxneEj=1…P求解以上P个方程组则可得到预测系数ak。方程组的求解方法有自相关的算法，有协方差的方法，本实验采用Durbin算法来求解，完整的递推过程如下：(1))0(0nnRE=(2)1111)()(−−=−∑−−=inikiknniEakiRiRk∑=−−=′−=pkkknxanxnxnxne1)()()()()((3)iiika=(4)iikiikikkaaa11−−−−=11−=ikL(5)12)1(−−=