clear;
x=wavread('a1.wav'); %读取声音文件
figure(1);stem(x,'.'); title('原始声音波形'); %显示声音波形
x=x(7701:13200); %采样频率11KHZ，共2.5s，共27500个采样点，截取其中0.7s到1.2s这一时间段，即从7701点到13200点，共5500个采样点;取帧长为20ms，即每帧220个采样点，共分为25帧。
figure(2);stem(x,'.'); title('截取的声音波形'); %显示声音波形
%运用短时自相关函数法检测基音周期
n=220; %取20ms的声音片段，即220个样点为一帧
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(3);stem(T,'.');axis([0 length(T) 0 20]);
xlabel('帧数(n)');ylabel('周期(ms)');title('各帧基音周期');
T1= medfilt1(T,5); %去除野点(基因周期偏离正常估计值的点)，使用中值平滑算法
figure(4);stem(T1,'.');axis([0 length(T1) 0 20]);
xlabel('帧数(n)');ylabel('周期(ms)');title('去除野点之后各帧基音周期');
