%%corr.m
% correlation for pitch estimation
% flag =1 from left to right, otherwise, from right to left
% r=corr(s,L,flag)
%
%
%flag==1时
%
%
%
%|--i----|
%s:.....................................................
% .
% .
% .
% .
% .
% .
%s:.....................................................
% |_______d________|__________________L-d_____________|
%返回值是一个序列
%注意到在做自相关时，错位值从20开始
%考虑到在用这个自相关序列时可能是用于求周期，这个20时有利的
%语音的基音周期一般在20到140，因此小于20没有周期
%而错开20以上可以避免出现不是所求的峰值避免出错
%求到结果后放到rr序列中时，注意到下标为i时实际上已经是错开
%i+20-1,因此pitch.m在调用这个时需要有一句p=19+ind_pitch;
function r=corr(s,d,flag)
epsilon=1.0e-10; % to avoid divided by zero
L=length(s);
NM=120;
if flag==1 % From left to right
for i=20:min(L-d,NM)
temp=0.0;
temp1=0.0;
temp2=0.0;
for j=1:d
temp=temp+s(j)*s(j+i);
temp1=temp1+s(j+i)*s(j+i);
temp2=temp2+s(j)*s(j);
end
pw1=sqrt(temp1)+epsilon;
pw2=sqrt(temp2)+epsilon;
r(i-19)=temp/(pw1*pw2);
end
else % From right to left
for i=20:min(L-d,NM)
temp=0.0;
temp1=0.0;
temp2=0.0;
for j=L-d+1:L
temp=temp+s(j-i)*s(j);
temp1=temp1+s(j-i)*s(j-i);
temp2=temp2+s(j)*s(j);
end
pw1=sqrt(temp1)+epsilon;
pw2=sqrt(temp2)+epsilon;
r(i-19)=temp/(pw1*pw2);
end
end
%%maxx.m
% Find all the local maxima of a function which are greater than
% t times the global maximum
% max_index=maxx(s,t)
%输入一个序列和一个标量t
%首先求出序列s的最大值max_value
%然后如果序列中的一个点比相邻两个点都大
%而且还大于max_value的t倍
%就把这个点的位置存放在序列mi中
%在序列mi最后再补上序列的最后一个点
%考虑到第一个点和最后一个点有可能满足要求
%但这两个点都不能跟旁边的比较
%所以还是把这两个点保留了，放在mi序列的第一个和最后一个
%实际上就是以最大值点的t倍为界找出这个界以上的极值点的位置
function mi=maxx(s,t)
max_value=max(s);
oldstate=0;
L=length(s);
j=1;
mi=[1];
for i=1:L-1
if s(i)>s(i+1)
newstate=-1;
end
if s(i)t*max_value
j=j+1;
mi(j)=i;
end
end
oldstate=newstate;
end
ll=length(mi);
mi(ll+1)=length(s);
%%pitch.m
% This is a function to find the pitch period of a long enough speech
% interval. In order to find correct pitch period, it is in favor of
% selecting short period. i.e., if p1p2*Thr, we select p1
% instead of p2.
%
% usage [p,mx]=pitch(s)
%输入一个序列s,首先求出其自相关序列rr
%要求出第一个周期的极大值点
%首先求出序列的最大值maxim所在的位置，如果这个位置在THR_largepitch以内
%就认为这就是要求的点。如果这个位置大于THR_largepitch，就考虑小于
%THR_largepitch的范围内有没有点满足以下条件：大于maxim*THR_MAXX的极值点，
%大于max(THR_corr,THR_pitch*maxim).
%如果有满足条件的，第一个点就是所求点
%至于为什么可以这样做，还没搞懂
%总的说来，这个函数返回一个序列的第一个周期
function [p,mx]=pitch(s)
THR_pitch=0.75;
THR_maxx=0.7;
THR_corr=0.33;
THR_largePitch=60;
epsilon=1.0e-10; % used to avoid dividing by zero
L=length(s);
d=100;
rr=corr(s,d,0);%最后一个参数是0，求自相关时从右到左
[maxim,mxi]=max(rr);
%当rr是一个行向量时，[Y,I]=max(s);Y是s中的最大值，I是标量
%存放了最大值的位置
ind_pitch=mxi;
if mxi>THR_largePitch
mi=maxx(rr,THR_maxx);%实际上就是以最大值点的t倍为界找出这个界以上的左右极值点的位置
%这个序列的第一个值是1，最后一个是序列rr的长度
%一开始看到这里的时候会觉得比较多余了，直接一次maxim*THR_pitch分界比较就可以了
%事实上，下面的比较没有只针对极值点处理。
%maxx的作用是筛选出大于maxim*THR_maxx的极值点出来，所以这个还是必须的
for i=1:length(mi)
if rr(mi(i))>max(THR_corr,maxim*THR_pitch);
ind_pitch=mi(i);
break;%break说明只取第一个满足条件的
end
end
end
p=19+ind_pitch;%这个19参考corr.m中的解释，其实就是rr(i)对应的是i+19
mx=rr(ind_pitch);
%%zerocros.m
% This is a function to calculate the zero-cross of a signal
% used to determine a signal is voiced or not.
%
% usage zc=zerocros(s)
%输入一个序列，如果相邻两个点异号，统计数字加1
%这个结果可以作为判断是否是噪音的一个依据
%当count足够大时就可以认为是噪音
function zc=zerocros(s)
L=length(s);
count=0;
for i=1:L-1
if s(i)*s(i+1)<0
count=count+1;
end
end
zc=count;
