multicuef0.m

该程序是日本人用matlab写出来的提取pitch的

smmp=fftfilt((hanning(sml).^2)/sum((hanning(sml).^2)),[mmp zeros(fftl/2+1,sml*2)]'+0.00001)';
smmp(smmp==0)=smmp(smmp==0)+0.0000000001;
%%waitbar(0.2);
smmp=1.0./fftfilt(hanning(sml)/sum(hanning(sml)),1.0./smmp')';
%%waitbar(0.4);
smmp=smmp(:,max(1,(1:nfr)+sml-2)); % fixed by H.K. on 10/Dec./2002

%--- Power adaptive weighting (29/July/1999)

spwm=fftfilt(hanning(sml)/sum(hanning(sml)),[pwm zeros(fftl/2+1,sml*2)]')';
spwm(spwm==0)=spwm(spwm==0)+0.00000001;
%%%waitbar(0.6);
spfm=fftfilt(hanning(sml)/sum(hanning(sml)),[pwm.*pif zeros(fftl/2+1,sml*2)]'+0.00001)';
%%%waitbar(0.8);
spif=spfm./spwm;
spif=spif(:,(1:nfr)+smb);

%keyboard;

idx=max(0,f0i/fs*fftl);
%iidx=idx+(0:nfr-1)*(fftl/2+1)+1;
%vv=smmp(floor(iidx))+(iidx-floor(iidx)).*(smmp(floor(iidx)+1)-smmp(floor(iidx)));
%iidx2=iidx+idx;
%iidx3=iidx+2*idx;
%vv2=smmp(floor(iidx2))+(iidx2-floor(iidx2)).*(smmp(floor(iidx2)+1)-smmp(floor(iidx2)));
%vv3=smmp(floor(iidx3))+(iidx3-floor(iidx3)).*(smmp(floor(iidx3)+1)-smmp(floor(iidx3)));
%fq1=(pif(floor(iidx))+(iidx-floor(iidx)).*(pif(floor(iidx)+1)-pif(floor(iidx))))/2/pi;
%fq2=(pif(floor(iidx2))+(iidx2-floor(iidx2)).*(pif(floor(iidx2)+1)-pif(floor(iidx2))))/2/pi;
%fq3=(pif(floor(iidx3))+(iidx3-floor(iidx3)).*(pif(floor(iidx3)+1)-pif(floor(iidx3))))/2/pi;
%ffq=(fq1./sqrt(vv)+fq2./sqrt(vv2/4)/2+fq3./sqrt(vv3/9)/3)./(1.0./sqrt(vv)+1.0./sqrt(vv2/4)+1.0./sqrt(vv3/9));
%vvt=1.0./(1.0./vv+1.0./vv2*4+1.0./vv3*9);

fqv=zeros(nhmx,nfr);
vvv=zeros(nhmx,nfr);

iidx=(0:nfr-1)*(fftl/2+1)+1;
for ii=1:nhmx
	iidx=idx+iidx;
	vvv(ii,:)=(smmp(floor(iidx))+(iidx-floor(iidx)).*(smmp(floor(iidx)+1)-smmp(floor(iidx))))/(ii*ii);
%	fqv(ii,:)=(pif(floor(iidx))+(iidx-floor(iidx)).*(pif(floor(iidx)+1)-pif(floor(iidx))))/2/pi/ii;
	fqv(ii,:)=(spif(floor(iidx))+(iidx-floor(iidx)).*(spif(floor(iidx)+1)-spif(floor(iidx))))/2/pi/ii; % 29/July/199
end;
%%%waitbar(1);
vvvf=1.0./sum(1.0./vvv);
f0r=sum(fqv./sqrt(vvv))./sum(1.0./sqrt(vvv)).*(f0raw>0);
ecr=sqrt(1.0./vvvf).*(f0raw>0)+(f0raw<=0);
%%%close(hpg);

%keyboard;

%--------------------
function [c1,c2]=znrmlcf3(f)

n=100;
x=0:1/n:3;
g=GcBs(x,0);
dg=[diff(g) 0]*n;
dgs=dg/2/pi/f;
xx=2*pi*f*x;
c1=sum((xx.*dgs).^2)/n;
c2=sum((xx.^2.*dgs).^2)/n;

%---------------------
function p=GcBs(x,k)

tt=x+0.0000001;
p=tt.^k.*exp(-pi*tt.^2).*(sin(pi*tt+0.0001)./(pi*tt+0.0001)).^2;

%%%----------------
function [f0,pl]=zcombineRanking4(p,beta,wAC,wIF,prm)

%        F0candidatesByIF: [2978x3 double]
%      CNofcandidatesByIF: [2978x3 double]
%        F0candidatesByAC: [2977x3 double]
%      ACofcandidatesByAC: [2977x3 double]
%               RefinedCN: [1x2977 double]
%    FirstAutoCorrelation: [1x2977 double]
%       F0initialEstimate: [2977x1 double]
%    BackgroundNoiselevel: -69.1041
%      InstantaneousPower: [1x2979 double]

f0floor=prm.F0searchLowerBound; % f0floor
f0ceil=prm.F0searchUpperBound; % f0ceil

n=min([length(p.F0candidatesByIF) length(p.F0candidatesByAC)]);
nr=n;

nvo=24;
nvc=ceil(log2(f0ceil/f0floor))*nvo;
fx=f0floor*2.0.^((0:nvc-1)/nvo);
lfx=log2(fx);
v1=[1;1;1];
r1=[1 1 1];
logf0if=log2(p.F0candidatesByIF);
logf0ac=log2(p.F0candidatesByAC);
relif=max(0.000000001,(p.CNofcandidatesByIF-min(p.CNofcandidatesByIF(:,1)))./ ...
    (max(p.CNofcandidatesByIF(:,1))-min(p.CNofcandidatesByIF(:,1))));
relac=max(0.000000001,((p.ACofcandidatesByAC-min(p.ACofcandidatesByAC(:,1)))./ ...
    (max(p.ACofcandidatesByAC(:,1))-min(p.ACofcandidatesByAC(:,1)))));
f0=zeros(n,6);
pl=zeros(n,6);
initv=0*lfx;
for ii=1:n
    IFmap=initv;
    ACmap=initv;
    for jj=1:3
        IFmap=IFmap+relif(ii,jj)^2*exp(-((logf0if(ii,jj)-lfx)/beta).^2); % 27/July/2004 ^2
        ACmap=ACmap+relac(ii,jj)^2*exp(-((logf0ac(ii,jj)-lfx)/beta).^2); % 27/July/2004 ^2
    end;
    %f0map=wIF*IFmap+wAC*ACmap;% 22/July/2004
    %f0map=sqrt(IFmap.*ACmap); % 22/July/2004 sqrt
    f0map=sqrt(wIF*IFmap+wAC*ACmap)/sqrt(2); % 27/July/2004 Addition
    f0mapbak=f0map;
    ix=find((diff([f0map(1) f0map]).*diff([f0map f0map(end)]))<0);
    if length(ix)>0
        [mx,mxp]=max(f0map(ix));
        [pl(ii,1),pos]=zzParabolicInterp2(f0mapbak((-1:1)+ix(mxp)),ix(mxp));
        [tpl,idsrt]=sort(-f0map(ix));
        nix=length(ix);
        for jj=1:6
            if jj>nix
                pl(ii,jj)=pl(ii,jj-1);
                f0(ii,jj)=f0(ii,jj-1);
            else
                [pl(ii,jj),f0(ii,jj)]=zzParabolicInterp2(f0mapbak((-1:1)+ix((idsrt(jj)))),ix((idsrt(jj))));
            end;
        end;
    else
        pos=round(nr/4);
        pl(ii,1)=0;
    end;
end;
f0=f0floor*2.0.^(f0/nvo);

%-------
function [val,pos]=zzParabolicInterp2(yv,xo)
lp=diff(yv);
a=lp(1)-lp(2);
b=(lp(1)+lp(2))/2;
xp=b/a+xo;
val=yv(2)+0.5*a*(b/a)^2+b*(b/a);
pos=xp-1;

%%%---
function ac1=zfirstac(x,fs,ix,wlms);

wl=round(fs*wlms/1000);
fftl=2.0^ceil(log2(wl));
xx=x(:);
idx=ix+(1:wl)-round(wl/2);
xt=xx(min(length(xx),max(1,idx)));
if sum(abs(xt))<1e-20  % bug fix 11/Jan./2005
    xt=xt+randn(size(xt));
end;
fw=abs(fft(xt.*hanning(wl),fftl)).^2;
fx=(0:fftl-1)/fftl*fs;
[dmy,ixx]=min(abs(fx-4000));
if ixx>fftl/2
    ixx=fftl/2;
end;
c=cos(2*pi*fx(1:ixx)/(fx(ixx)*2));
ac1=sum(c'.*fw(1:ixx))/sum(fw(1:ixx));

%%%-------

function [f0raw,segb]=zf0trackmulti25(p,ac1th,f0jump,nsd,prm)

%       F0candidatesByMix: [2789x3 double]
%    RELofcandidatesByMix: [2789x3 double]
%               RefinedCN: [1x2789 double]
%    FirstAutoCorrelation: [1x2789 double]
%       F0initialEstimate: [2789x1 double]
%    BackgroundNoiselevel: -51.0845
%      InstantaneousPower: [1x2791 double]

f0floor=prm.F0searchLowerBound; % f0floor
f0ceil=prm.F0searchUpperBound; % f0ceil

pwsdb=p.InstantaneousPower;
f0cand=p.F0candidatesByMix;
relv=p.RELofcandidatesByMix;
relv(relv==0)=relv(relv==0)+0.00001;
acv=p.FirstAutoCorrelation;
nn=min(length(pwsdb),length(f0cand));
DispOn=prm.DisplayPlots;
f0jumpt=f0jump;
nsdt=nsd;

if DispOn==1
    figure;
    semilogy(f0cand,'+');grid on;
    axis([0 length(f0cand) f0floor f0ceil]);
    drawnow
    hold on
end;

%---- noiselevel
%pwsdb=10*log10(pws+0.00000000001);
pws=10.0.^(pwsdb/10);
mxpwsdb=max(pwsdb);
[hstgrm,binlvl]=hist(pwsdb,mxpwsdb+(-60:2));
q10=interp1(cumsum(hstgrm+0.000000001)/sum(hstgrm)*100,binlvl,10); % 10% quantile level
[dmy1,minid]=min(abs(q10-binlvl));
%bb=minid+(-5:5); % search range 10 dB
bb=max(1,min(length(binlvl),minid+(-5:5))); % search range 10 dB % safeguard
noiselevel=sum(hstgrm(bb).*binlvl(bb))/sum(hstgrm(bb)); 

%----- mark syllable centers
wsml=81;
pwsdbl=[ones(wsml,1)*pwsdb(1);pwsdb(:);ones(2*wsml,1)*pwsdb(end)];
pwsdbs=fftfilt(hanning(wsml)/sum(hanning(wsml)),pwsdbl);
pwsdbs=pwsdbs((1:length(pwsdb))+round(3*wsml/2));
dpwsdbs=diff([pwsdbs(1);pwsdbs]);
dpwsdbsm=diff([pwsdbs;pwsdbs(end)]);
pv=find((dpwsdbs.*dpwsdbsm<0)&(dpwsdbsm<=0));
dv=find((dpwsdbs.*dpwsdbsm<0)&(dpwsdbsm>0));
%pv=pv(pwsdbs(pv)>(noiselevel+3));
% -------

[nr,nc]=size(f0cand);
nn=min(nr,length(pws));
dv=[1;dv(:);nn]; % safe guard
f0raw0=zeros(nn,1);
relf0=f0raw0;
%f0jump=0.3;
%kk=0;
%lastsegf0=0;
%for ii=1:length(pv)
%    ac1=zfirstac(x,fs,round(pv(ii)/1000*fs),30);
%    if ac1>0.4
%        kk=kk+1;
%        lastsegf0=f0cand(pv(ii),1)+lastsegf0;
%    end;
%end;
%lastsegf0=lastsegf0/kk
%lastsegf0=sum(pws(1:nn)'.*f0cand(1:nn,1))/sum(pws(1:nn));

lastsegf0=2.0^(sum(pws(1:nn)'.*log2(f0cand(1:nn,1)))/sum(pws(1:nn)));

lastpos=1;
segb=zeros(length(pv),1);
finalf0=lastsegf0;
for ii=1:length(pv)
    %    [mxv,mxp]=max(relv((-10:10)+pv(ii),1)); % This OPs may be fragile.
    %    fqs=sort(f0cand((-10:10)+pv(ii),1));
    %    [mxv,mxp]=min(abs(fqs(10)-f0cand((-10:10)+pv(ii),1)));
    %    acp=mxp+pv(ii)-11;
    acp=pv(ii);
    lb=max(dv(dv<acp));
    if lb>lastpos
        lb=lastpos;
    end;
    ub=min(dv(dv>acp));
    if length(lb)==0;lb=1;end;
    if length(ub)==0;ub=nn;end;
    fqs=sort(f0cand(lb:ub,1));
    medianf0=fqs(round((ub-lb)/2));
    [mxv,mxp]=min(abs(lastsegf0-f0cand(max(lb,min(ub,(-10:10)+pv(ii))),1)));
    acp=max(lb,min(ub,mxp+pv(ii)-11));
    f0raw0(acp)=f0cand(acp,1);%f0cand(acp,1);
    lastf0=f0raw0(acp);
    %    ac1=zfirstac(x,fs,round(pv(ii)/1000*fs),30);
    ac1=acv(acp);
    %    lb
    %    if (ub-lb)>12; 
    %        std(diff(log2(f0cand((lb+5):(ub-5),1))))
    %    end;
    segb(ii)=lb;
    forwardind=1;
    if DispOn==1
        semilogy([lb lb],[f0floor f0ceil]);
        semilogy([ub ub],[f0floor f0ceil],'r');
        hf=semilogy(acp,lastf0,'ro');
        set(hf,'MarkerSize',10,'LineWidth',2)
    end;
    %[pv(ii),lb,ub,acp]
    %    lastsegf0=lastf0;
    f0bak=f0raw0;
    f0raw1=f0raw0;
    f0raw2=f0raw0;
    f0raw3=f0raw0;
    [f0rawm,sprob0]=ztraceInAsegment(f0raw0,f0cand,relv,pwsdb,acp,lastf0,lb,ub,nn,f0jumpt,nsdt,noiselevel);
    f0raw1(acp)=f0cand(acp,2);%f0cand(acp,1);
    lastf0=f0raw1(acp);
    [f0raw1,sprob1]=ztraceInAsegment(f0raw1,f0cand,relv,pwsdb,acp,lastf0,lb,ub,nn,f0jumpt,nsdt,noiselevel);
    f0raw2(lb)=f0cand(lb,1);%f0cand(acp,1);
    lastf0=f0raw2(lb);
    [f0raw2,sprob2]=ztraceInAsegment(f0raw2,f0cand,relv,pwsdb,lb,lastf0,lb,ub,nn,f0jumpt,nsdt,noiselevel);
    f0raw2(ub)=f0cand(ub,1);%f0cand(acp,1);
    lastf0=f0raw2(ub);
    [f0raw3,sprob3]=ztraceInAsegment(f0raw3,f0cand,relv,pwsdb,ub,lastf0,lb,ub,nn,f0jumpt,nsdt,noiselevel);
%    f0raw0=f0bak;
%    sprob1=0;
    %[[lb ub]/1000 sprob0 sprob1 sprob2 sprob3]
    [mx,imx]=max([sprob0 sprob1 sprob2 sprob3]);
%    if sprob1>sprob0
%        f0raw0=f0raw1;
%    end;
    switch imx
        case 1
            f0raw0=f0rawm;
        case 2
            f0raw0=f0raw1;
        case 3
            f0raw0=f0raw2;
        case 4
            f0raw0=f0raw3;
    end;
    %[lb ub lastpos finalf0  f0raw0(lb)]
    reliablepowerth=max([noiselevel+10, min([(2*mxpwsdb+noiselevel)/3])]);%, mxpwsdb+max(pwsdb(lb:ub))-5])]); %24/July/2004 #2
    %[lb ub mxpwsdb noiselevel max(pwsdb(lb:ub)) reliablepowerth]
    if (abs(log2(finalf0)-log2(f0raw0(lb)))>0.4)&(ii>1)&(pwsdb(lb)>reliablepowerth) %&(0==1) % (mxpwsdb+noiselevel)/2)
        disp(['hit! at ' num2str(lb)]);
        f0raw4=f0raw0;
        lastf0=f0raw0(lb);
        [f0raw4,sprob4]=ztraceInAsegment(f0raw4,f0cand,relv,pwsdb,lb,lastf0,segb(ii-1),lb,nn,f0jumpt,nsdt,noiselevel);
        %sprob4
        %[ finalf0 f0rawm(lb) f0raw1(lb) f0raw2(lb) f0raw3(lb)]
        %abs(log2(finalf0)-log2([f0rawm(lb) f0raw1(lb) f0raw2(lb) f0raw3(lb)]))
        [secondCand,idsec]=min(abs(log2(finalf0)-log2([f0rawm(lb) f0raw1(lb) f0raw2(lb) f0raw3(lb)])));
        switch idsec
            case 1
                secprob=sprob0;
                f0sec=f0rawm;
            case 2
                secprob=sprob1;
                f0sec=f0raw1;
            case 3
                secprob=sprob2;
                f0sec=f0raw2;
            case 4
                secprob=sprob3;
                f0sec=f0raw3;
        end;
        jprob1=exp(-((log2(finalf0)-log2(f0raw0(lb)))/nsd)^2);
        jprob2=exp(-(secondCand/nsd)^2);
        %[mx*lastprob*jprob1 sprob4*mx lastprob*jprob2*secprob]
        [mx2,imx2]=max([mx*lastprob*jprob1 sprob4*mx lastprob*jprob2*secprob]);
        switch imx2
            case 1
                % do nothing
            case 2
                f0raw0(segb(ii-1):lb)=f0raw4(segb(ii-1):lb);
            case 3
                f0raw0(lb:ub)=f0sec(lb:ub);
        end;
    end;
    lastpos=ub;
    lastprob=mx;
    finalf0=f0raw0(ub);
end;
segb=[segb;nn];
f0raw=f0raw0;

%----
function [f0raw0,sprob]=ztraceInAsegment(f0rawin,f0cand,relv,pwsdb,acp,lastf0in,lb,ub,nn,f0jump,nsd,noiselevel)

forwardind=1;
lastf0=lastf0in;
f0raw0=f0rawin;
maxpower=max(pwsdb); %(pwsdb(acp)+max(pwsdb))/2;
%reliablepowerth=maxpower-12;
%reliablepowerth=(maxpower+3*noiselevel)/4;
reliablepowerth=(2*maxpower+noiselevel)/3; % 23/July/2004
reliablepowerth=(4*pwsdb(acp)+noiselevel)/5; % 24/July/2004
reliablepowerth=max([noiselevel+10, min([(3*maxpower+noiselevel)/4])]);%, max(pwsdb(lb:ub))-5])]); %24/July/2004 #2
sprob=0; % initial probability
if acp==ub;f0raw0(ub)=lastf0in;end;
if acp==lb;f0raw0(lb)=lastf0in;end;
for jj=acp-1:-1:lb
    bsb=max(1,jj:-1:jj-3); % seraching bound
    [dmy,idx]=max(exp(-((log2(f0cand(bsb,:)')-log2(lastf0))/nsd).^2).*(relv(bsb,:)'));
    [dmy2,idxx]=max(dmy);
    idx=idx(idxx);
    jjmx=bsb(idxx);
    if abs(log2(f0cand(jjmx,idx))-log2(lastf0))<f0jump
        dd=abs(jjmx-jj);
        if dd==0;
            f0raw0(jj)=f0cand(jjmx,idx);
            if pwsdb(jj)>reliablepowerth;sprob=sprob+log2(dmy2);end;
        else 
            f0raw0(jj)=lastf0+(1/(dd+1))*(f0cand(jjmx,idx)-lastf0);
            if pwsdb(jj)>reliablepowerth;sprob=sprob+log2(dmy2);end;%-0.1;end;
        end;
    else
        f0raw0(jj)=lastf0;
        if pwsdb(jj)>reliablepowerth;sprob=sprob+log2(dmy2)-10;end;
    end
    lastf0=f0raw0(jj);
end;
lastf0=lastf0in;%f0cand(acp,1);
for jj=acp+1:ub
    bsb=min(nn,jj:jj+3); % seraching bound
    [dmy,idx]=max(exp(-((log2(f0cand(bsb,:)')-log2(lastf0))/nsd).^2).*(relv(bsb,:)'));
    [dmy2,idxx]=max(dmy);
    idx=idx(idxx);
    jjmx=bsb(idxx);
    if abs(log2(f0cand(jjmx,idx))-log2(lastf0))<f0jump
        dd=abs(jjmx-jj);
        if dd==0;
            f0raw0(jj)=f0cand(jjmx,idx);
            if pwsdb(jj)>reliablepowerth;sprob=sprob+log2(dmy2);end;
        else 
            f0raw0(jj)=lastf0+(1/