www.gusucode.com > matlab编程倒谱法来进行基音周期检测源码程序 > gusucode/CEP/CEP2.m

      clear all
  load ma1_4;
  x=ma1_4;
  figure,plot(x);
  len=length(x);
  fss=enframe(x,256,128);xlabel('采样点数');ylabel('幅度');title('原始语音图像');
  %分帧加窗
  F=8000;
  [row,col]=size(fss);%记录分帧后矩阵大小：row代表帧数，col代表帧长
  T=zeros(1,row);%用来存储检测出的基音周期对应的抽样点数
  zeronum=0;sum1=0;k=0;
  zcr=countoverzeroing(x,256,128);      
  %计算过零率
  figure,plot(zcr,'.');xlabel('分帧数');ylabel('过零率');title('过零率统计');grid on;
  %绘制过零率图像
  for m=1:row-1
  f1=zeros(1,col);%用来存储当前正在检测基音周期的帧
  for n=1:col
   f1(1,n)=fss(m,n);
  end
  A1=lpc(f1,18);
  E1=abs(filter(A1,1,f1));
% LPC系数 逆滤波
  e1=abs(fftshift(fft(E1)));
% FFT
if e1~=0
  e1=log(e1);
end
% 绝对值取对数 
   for n=1:128
       e1(n)=0;
   end
   for n=144:256
        e1(n)=0;
   end
%高频分量置零 500/800*128+128=144
 EE1=abs(ifft(e1));
%  IDFT
 for n=1:64
     EE1(n)=0;
 end
 for n=181:256
     EE1(n)=0;
 end
%噪音尖峰置零
  N=max(EE1);
 for n=1:128
     if abs(EE1(n)-N)<0.0001
         n1=n;
     end
 end
%根据峰值，检出基音周期
  if zcr(m,1)>1 && zcr(m,1)<8
      T(1,m)=n1;sum1=sum1+T(1,m);k=k+1;
  else
       T(1,m)=0;zeronum=zeronum+1;
  end  
  end
%根据过零率，去除静音，清音
mean=sum1/(row-zeronum); 
zeronum=0;sum=0;
for m=1:row
    if T(1,m)<mean*0.5|T(1,m)>mean*1.5
               T(1,m)=0;zeronum=zeronum+1;
    end
    sum=sum+T(1,m);
end   
turemean=sum/(row-zeronum);
%去除幅度与均值相差较大的野点,剔除后再求平均值
 zeronum=0;sum=0;
 for m=1:row
    if T(1,m)<turemean*0.8|T(1,m)>turemean*1.2
               T(1,m)=F;zeronum=zeronum+1;
    end
    sum=sum+T(1,m);
    if T(1,m)~=0
     T(1,m)=1/T(1,m);
    end
end
%第二次去除野点,求平均值
turemean1=(sum-zeronum*F)/(row-zeronum);
figure,plot(F*T,'.');hold on;xlabel('分帧数');ylabel('基音频率/HZ');title('基音周期统计');
%转换各帧周期为频域数据,绘图
for m=1:row
  if turemean1~=0    
    T(1,m)=1/turemean1; 
  end
end  
  plot(F*T,'r');hold off;
  axis([0 row+5 0 140]);grid on;
%绘制基音周期平均值
   fmean=F*T(1)