语音基音周期的估计.doc

实验一 离散数字信号的产生及其时域处理 学习实现实验1的内容，并且编制一个程序(m文件)产生5种信号，函数需要的参数可输入确定，并绘出其图形1、单位抽样序列 在MATLAB中可以利用函数实现,%单位抽样序列函数% X=0:10;Y=0 1 zeros(1,9);stem(X,Y,r);axis(-1,10,0,1);title(单位抽样序列 );xlabel(n);ylabel(n);图形如右: 2、单位阶越序列 在MATLAB中可以利用函数实现, 实现过程如下：%单位阶跃序列函数K=-8:8;H=zeros(1,8),ones(1,9);stem(K,H,r);axis(-8,8,0,2);title(单位阶跃序列 );xlabel(n);ylabel(un);图形如下:3、正弦序列 ，在MATLAB中实现过程如下：%正弦序列函数sin(2*pi*D/5+pi/4)%D=-1:0.1*pi:8*pi;C=sin(2*pi*D/5+pi/4);stem(D,C,filled);axis(-1,10,-2,2);title(正弦序列 );xlabel(n);ylabel(sin(2*pi*D/5+pi/4)图形如下:4、复指数序列 ，从幅度和相位进行分析，在MATLAB中实现过程如下：%复指数序列函数%n=0:10;x1=2*exp(-0.2+0.7*j)*n);x2=abs(x1);x3=angle(x1);subplot(2,3,4);stem(n,x2);title(复指数序列幅值);xlabel(n);ylabel(x2);stem(n,x3);title(复指数序列相位);xlabel(n);ylabel(x3);图形如下：5、指数序列 ，在MATLAB中实现过程如下： % 实指数序列%L=0:20;x4=(1.4).(L/3);stem(L,x4);title(实指数序列 );xlabel(n);ylabel(x4);图形如下：编制一个完整的程序(m文件)产生以上5种信号，图形如下：m文件程序如下：function F1%单位抽样序列函数% X=0:10;Y=0 1 zeros(1,9);subplot(2,3,1); stem(X,Y,r);axis(-1,10,0,1);title(单位抽样序列 );xlabel(n);ylabel(n);%单位阶跃序列函数K=-8:8;H=zeros(1,8),ones(1,9);subplot(2,3,2);stem(K,H,r);axis(-8,8,0,2);title(单位阶跃序列 );xlabel(n);ylabel(un);%正弦序列函数sin(2*pi*D/5+pi/4)%D=-1:0.1*pi:8*pi;C=sin(2*pi*D/5+pi/4);subplot(2,3,3);stem(D,C,filled);axis(-1,10,-2,2);title(正弦序列 );xlabel(n);ylabel(sin(2*pi*D/5+pi/4);%复指数序列函数2*exp(-0.2+0.7*j)*n)，从相位和幅值角度来分开讨论其图形%n=0:10;x1=2*exp(-0.2+0.7*j)*n);x2=abs(x1);x3=angle(x1);subplot(2,3,4);stem(n,x2);title(复指数序列幅值);xlabel(n);ylabel(x2);subplot(2,3,5);stem(n,x3);title(复指数序列相位);xlabel(n);ylabel(x3);% 实指数序列(1.4)(L/3)%L=0:20;x4=(1.4).(L/3);subplot(2,3,6);stem(L,x4);title(实指数序列 );xlabel(n);ylabel(x4);%在MATLAB中产生5种信号，所需参数预先设定% 语音基音周期估计 语音基音周期估计的实现方法: 自相关函数法能量有限的语音信号的短时自相关函数定义为： 其中，为移位距离，是偶对称的窗函数。短时自相关函数有以下重要性质：如果是周期信号，周期是，则也是周期信号，且周期相同，即。当=0时，自相关函数具有最大值；当处周期信号的自相关函数达到极大值。自相关函数是偶函数，即。短时自相关函数法基音检测的主要原理是利用短时自相关函数的第二条性质，通过比较原始信号和它移位后的信号之间的类似性来确定基音周期，如果移位距离等于基音周期，那么，两个信号具有最大类似性。在实际采用短时自相关函数法进行基音检测时，使用一个窗函数，窗不动，语音信号移动，这是经典的短时自相关函数法。窗口长度的选择至少要大于基音周期的两倍，越大，短时自相关函数波形的细节就越清楚，更有利于基音检测，但计算量较大，近年来由于高速数字信号处理器(DSP)的使用，从而使得这一算法简单有效，而不再采用结构复杂的快速傅里叶变换法、递归计算法等；越小，误差越大，但计算量较小。自相关函数在基音周期处表现为峰值，自相关函数在基音周期处表现为峰值，这些峰值点之间的间隔的平均值就是基音周期中值滤波: 为了平滑噪声，君安邪恶基因周期前，需要对各自的声音文件进行中值滤波处理，并比较前后语音的差别。对采样信号进行滤波处理，这里用截至频率为500Hz的低通滤波器,阻带衰减20dB，程序如下：x,fs,bits=wavread(E:A.wav);%读取输入语音信号的频率及波形figure(1);stem(x,.);%做原始语音信号的时域图形title(原始语音信号);fc1=500;N1=2*pi*0.9/(0.1*pi)wc1=2*pi*fc1/fs;if rem(N1,2)=0 N1=N1+1;endWindow= boxcar (N1+1); %长度为N1的矩形窗Windowb1=fir1(N1,wc1/pi,Window);figure(2);freqz(b1,1,512);title(低通滤波器的频率响应);y= filter(b1,1,x1);%对信号进行低通滤波figure(3);plot(y);title(信号经过低通滤波器后); 经过中值滤波后，对信号用自相关函数法进行基音周期的检测。语音信号需分帧处理，20ms一帧（160个样点） ，程序如下： n=160; %取20ms的声音片段，即160个样点for m=1:length(x)/n; %对每一帧求短时自相关函数 for k=1:n; Rm(k)=0; for i=(k+1):n; Rm(k)=Rm(k)+x(i+(m-1)*n)*x(i-k+(m-1)*n); end end p=Rm(10:n); %防止误判，去掉前边10个数值较大的点 Rmax,N(m)=max(p); %读取第一个自相关函数的最大点end %补回前边去掉的10个点N=N+10;T=N/8; %算出对应的周期 figure(4);stem(T,.);axis(0 length(T) 0 10);xlabel(帧数(n);ylabel(周期(ms);title(各帧基音周期);得到各帧基音周期如图：可以看到， 该信号有48000个样点，一共有300帧（每帧160个样点）由图中可以看出基音周期大约为4.2ms.但是图中存在野点，为此，需要对此进行进一步的处理，即去除野点。运行以下程序：T1= medfilt1(T,5); %去除野点figure(3);stem(T1,.);axis(0 length(T1) 0 20); xlabel(帧数(n);ylabel(周期(ms);title(各帧基音周期);从图可以看到野点被消除，可以读出基音周期约为4.2ms.完整程序如下：x,fs,bits=wavread(E:A.wav);figure(1);plot(x);%做原始语音信号的时域图形title(原始语音信号);fc1=500;N1=2*pi*0.9/(0.1*pi)wc1=2*pi*fc1/fs;if rem(N1,2)=0 N1=N1+1;endWindow= boxcar (N1+1); %长度为N1的矩形窗Windowb1=fir1(N1,wc1/pi,Window);figure(2);freqz(b1,1,512);title(低通滤波器的频率响应);y= filter(b1,1,x1);%对信号进行低通滤波figure(3);plot(y);title(信号经过低通滤波器(时域);n=160; %取20ms的声音片段，即160个样点for m=1:length(x)/n; %对每一帧求短时自相关函数 for k=1:n; Rm(k)=0; for i=(k+1):n; Rm(k)=Rm(k)+x(i+(m-1)*n)*x(i-k+(m-1)*n); end end p=Rm(10:n); %防止误判，去掉前边10个数值较大的点 Rmax,N(m)=max(p); %读取第一个自相关函数的最大点end %补回前边去掉的10个点N=N+10;T=N/8; %算出对应的