gusucode.com > matlab编程倒谱法来进行基音周期检测源码程序 > gusucode/CEP/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)