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)