cep1.m

倒谱法检测基音周期

clear all
load vowels;
x=vowels.a_2;
len=length(x);
fss=enframe(x,256,128);
%分帧加窗
F=10000;
[row,col]=size(fss);%记录分帧后矩阵大小：row代表帧数，col代表帧长
T=zeros(1,row);%用来存储检测出的基音周期对应的抽样点数
for m=1:row 
  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
  e1=log(e1);
% 绝对值取对数    
   for n=1:128
       e1(n)=0;
    end
    for n=139:256
        e1(n)=0;
    end
%高频分量置零 500/5000*128+128=140
 EE1=abs(ifft(e1));
 %  IDFT
 for n=1:79
       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
     T(1,m)=n1;    
 end
%根据峰值，检出基音周期
zeronum=0;sum=0;
 for m=1:row
    if T(1,m)<mean(T)*0.5|T(1,m)>mean(T)*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);
    T(1,m)=1/T(1,m);
end   
turemean1=(sum-zeronum*F)/(row-zeronum)
%第二次去除野点,求平均值
figure,plot(F*T,'.');hold on;xlabel('分帧数');ylabel('基音频率/HZ');title('基音周期统计');grid on;
%转换各帧周期为频域数据,绘图
for m=1:row
  T(1,m)=1/turemean1;    
end
plot(F*T,'r');hold off;
axis([0 row+5 0 250]);
%绘制基音周期平均值
fmean=F*T(1)