语音中提取基音频率matlab程序

%%corr.m%语音分帧，加窗%lile%correlationforpitchestimation%flag=1fromlefttoright,otherwise,fromrighttoleft%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;functionr=corr(s,d,flag)epsilon=1.0e-10;%toavoiddividedbyzeroL=length(s);NM=120;ifflag==1%Fromlefttorightfori=20:min(L-d,NM)temp=0.0;temp1=0.0;temp2=0.0;forj=1:dtemp=temp+s(j)*s(j+i);temp1=temp1+s(j+i)*s(j+i);temp2=temp2+s(j)*s(j);endpw1=sqrt(temp1)+epsilon;pw2=sqrt(temp2)+epsilon;r(i-19)=temp/(pw1*pw2);endelse%Fromrighttoleftfori=20:min(L-d,NM)temp=0.0;temp1=0.0;temp2=0.0;forj=L-d+1:Ltemp=temp+s(j-i)*s(j);temp1=temp1+s(j-i)*s(j-i);temp2=temp2+s(j)*s(j);endpw1=sqrt(temp1)+epsilon;pw2=sqrt(temp2)+epsilon;r(i-19)=temp/(pw1*pw2);endend%%maxx.m%Findallthelocalmaximaofafunctionwhicharegreaterthan%ttimestheglobalmaximum%max_index=maxx(s,t)%输入一个序列和一个标量t%首先求出序列s的最大值max_value%然后如果序列中的一个点比相邻两个点都大%而且还大于max_value的t倍%就把这个点的位置存放在序列mi中%在序列mi最后再补上序列的最后一个点%考虑到第一个点和最后一个点有可能满足要求%但这两个点都不能跟旁边的比较%所以还是把这两个点保留了，放在mi序列的第一个和最后一个%实际上就是以最大值点的t倍为界找出这个界以上的极值点的位置functionmi=maxx(s,t)max_value=max(s);oldstate=0;L=length(s);j=1;mi=[1];fori=1:L-1ifs(i)s(i+1)newstate=-1;endifs(i)t*max_valuej=j+1;mi(j)=i;endendoldstate=newstate;endll=length(mi);mi(ll+1)=length(s);%%pitch.m%Thisisafunctiontofindthepitchperiodofalongenoughspeech%interval.Inordertofindcorrectpitchperiod,itisinfavorof%selectingshortperiod.i.e.,ifp1p2*Thr,weselectp1%insteadofp2.%%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;%usedtoavoiddividingbyzeroL=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;ifmxiTHR_largePitchmi=maxx(rr,THR_maxx);%实际上就是以最大值点的t倍为界找出这个界以上的左右极值点的位置%这个序列的第一个值是1，最后一个是序列rr的长度%一开始看到这里的时候会觉得比较多余了，直接一次maxim*THR_pitch分界比较就可以了%事实上，下面的比较没有只针对极值点处理。%maxx的作用是筛选出大于maxim*THR_maxx的极值点出来，所以这个还是必须的fori=1:length(mi)ifrr(mi(i))max(THR_corr,maxim*THR_pitch);ind_pitch=mi(i);break;%break说明只取第一个满足条件的endendendp=19+ind_pitch;%这个19参考corr.m中的解释，其实就是rr(i)对应的是i+19mx=rr(ind_pitch);%%zerocros.m%Thisisafunctiontocalculatethezero-crossofasignal%usedtodetermineasignalisvoicedornot.%%usagezc=zerocros(s)%输入一个序列，如果相邻两个点异号，统计数字加1%这个结果可以作为判断是否是噪音的一个依据%当count足够大时就可以认为是噪音functionzc=zerocros(s)L=length(s);count=0;fori=1:L-1ifs(i)*s(i+1)0count=count+1;endendzc=count;