multicuef0.m

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

%   16/June/2004 Simplified version
%   17/June/2004 Spectral normalization version

%[x,fs]=wavread('m01.wav');
x=x(:);  % make x a column vector

wtlms=round(wtlm/1000*fs);  % windowlength in samples
wtlmso=floor(wtlms/2)*2+1;
bb=(1:wtlmso)-(wtlmso-1)/2; % time base for window;

fftl=2^ceil(log2(wtlmso)); % set FFT length to 2's exponent
x=[zeros(fftl,1);x;zeros(fftl,1)]; % safeguard

p=round(pm/1000*fs); % analysis position in samples
fx=(0:fftl-1)/fftl*fs;
tx=(0:fftl-1);
tx(tx>fftl/2)=tx(tx>fftl/2)-fftl;
tx=tx/fs;
lagw=exp(-(tx/0.004).^2);
lagw=exp(-(tx/0.003).^2); % EGGF0testn11
lagw=exp(-(tx/0.0035).^2); % EGGF0testn12

lagw2=exp(-(tx/0.0012).^2);
lagw2=exp(-(tx/0.0016).^2);% EGGF0testn12

xt=x(fftl+bb+p); % waveform segment to be analyzed
if sum(abs(xt))<1e-10 % bug fix 11/Jan./2005
    xt=xt+randn(size(xt));
end;
abase=abs(fft(xt.*blackman(wtlmso),fftl));
mxab=max(abase);
ac=ifft(abase.^2);
npw=real(fft(ac.*lagw'));
%efeps=min(npw(npw>0)); % bug fix 11/Jan./2005
%npw(npw==0)=npw(npw==0)+efeps; % bug fix 11/Jan./2005
pw=abase.^2.0./real(npw);

fsp=fs/fftl;
wflfs=round(wflf/fsp); % frequency window length in bins
wflfso=floor(wflfs/2)*2+1;
bbf=(1:wflfso)-(wflfso-1)/2; % index for frequency window
fftlf=2^ceil(log2(wflfso)+2);
lx=(0:fftlf/2-1)/(fsp*fftlf);
nsht=fftl/2/ndiv;
acc=zeros(fftlf/2,ndiv+1);
w2=hanning(wflfso);
ampw=1-lagw*(1-1/amp);
ampw=(1-lagw2(1:fftlf/2)'*(1-1/amp))./ampw(1:fftlf/2)';
for ii=1:ndiv+1
    p=rem(round(fftl/2+bbf+(ii-1)*nsht),fftl)+1;
    ac=abs(fft((pw(p)).*w2,fftlf))*(npw(p((wflfso-1)/2))).^pc;
    acc(:,ii)=ac(1:fftlf/2).*ampw;
end;

%%%----
function [f0,pl]=zmultiCand(lx,lagspec,beta,lagsp,timesp)
%   F0 candidates extraction from time-lag representation
%   using palabolic interpolation
%   [f0,pl]=multiCand(lx,lagspec,beta,lagsp,timesp)
%   lx  : lag axis
%   lagspec : time-lag representation
%   beta    : nonlinear distance measure
%   lagsp   : lag smoothing parameter (ms)
%   timesp  : temporal smoothing parameter (ms)
%   output parameters
%   f0  : fundamental frequency (Hz)
%   pl  : peak level

%   Designed and coded by Hideki Kawahara
%   30/March/2004
%   16/June/2004 Peak picking first
%   17/June/2004 Peak selection taking into interaction account

[nr,nc]=size(lagspec);
imm=diff([lagspec(1,:);lagspec]).*(diff([lagspec;lagspec(end,:)]));
dlag=diff([lagspec(1,:);lagspec]);
lagspecz=lagspec;
lagspecz(lx<0.001,:)=lagspec(lx<0.001,:)*0.0001;
%keyboard;
tls=[lagspecz;lagspecz(end,:);lagspecz(end:-1:2,:)].^beta;

llx=[lx lx(end) lx(end:-1:2)]';
lagw=exp(-(llx/(lagsp/1000)).^2); % This should be propotional to lag.
lagw=lagw/sum(lagw);
flagw=real(fft(lagw));
for ii=1:nc
    tls(:,ii)=real(ifft(fft(tls(:,ii)).*flagw));
end;

tmsm=round((timesp-1)/2)*2+1; % temporal smoothing (ms) assuming shift=1 ms
wt=hanning(tmsm);
wt=wt/sum(wt);
% temporal smoothing using 
lagsms=fftfilt(wt,[zeros(nr,tmsm) tls(1:nr,:) zeros(nr,tmsm)]')';
lagsms=lagsms(:,(1:nc)+(tmsm-1)/2*3);

lagsms=abs(lagsms).^(1/beta);
%figure;imagesc(lagsms);axis('xy');
f0=zeros(nc,3);
pl=zeros(nc,3);
for ii=1:nc
    ix=find((imm(:,ii)<0)&(dlag(:,ii)>0));
    [mx,mxp]=max(lagsms(ix,ii));
    [pl(ii,1),pos]=zParabolicInterp(lagspec((-1:1)+ix(mxp),ii),ix(mxp));
    f0(ii,1)=pos/lx(2);
    if length(ix)>1
        lagsms(ix(mxp),ii)=lagsms(ix(mxp),ii)*0;
        [mx,mxp]=max(lagsms(ix,ii));
        [pl(ii,2),pos]=zParabolicInterp(lagspec((-1:1)+ix(mxp),ii),ix(mxp));
        f0(ii,2)=pos/lx(2);
        if length(ix)>2
            lagsms(ix(mxp),ii)=lagsms(ix(mxp),ii)*0;
            [mx,mxp]=max(lagsms(ix,ii));
            [pl(ii,3),pos]=zParabolicInterp(lagspec((-1:1)+ix(mxp),ii),ix(mxp));
            f0(ii,3)=pos/lx(2);
        else
            f0(ii,3)=f0(ii,2);pl(ii,3)=pl(ii,2);
        end;
    else
        f0(ii,2)=f0(ii,1);pl(ii,2)=pl(ii,1);
    end;
end;

function [val,pos]=zParabolicInterp(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=1/(xp-1);

%%%-----
function [f0mapIF,f0mapAC,fx]=zCombineSourceInfoNorm(val,pos,f02,pl2,f0floor,f0ceil,nvo);
%   [f0mapIF,f0mapAC,fx]=CombineSourceInfo(val,pos,f02,pl2,f0floor,f0ceil,nvo);
%   val : C/N based on reliablity (dB)
%   pos : candidate frequency (Hz) based on Instantaneous frequency
%   f02 : Summary autocorrelation based candidate frequency (Hz)

[nrowIF,nc1]=size(val);
[nrowAC,nc2]=size(f02);

nn=min(nrowIF,nrowAC);
%f0floor=40;
%f0ceil=800;
%nvo=24;
nvc=ceil(log2(f0ceil/f0floor)*nvo);
fx=f0floor*2.0.^((0:nvc)/nvo)';
f0mapIF=zeros(length(fx),nn);
f0mapAC=zeros(length(fx),nn);

lb1=min(val(:,1));
lb2=min(pl2(:,1));
ub1=max(val(:,1));
ub2=max(pl2(:,1));

for ii=1:nn
    for jj=1:nc1
        f0mapIF(:,ii)=f0mapIF(:,ii)+((max(lb1,val(ii,jj))-lb1)/(ub1-lb1)).^2*exp(-(log2(pos(ii,jj)./fx)*6).^2);
        f0mapAC(:,ii)=f0mapAC(:,ii)+((max(lb2,pl2(ii,jj))-lb2)/(ub2-lb2)).^2*exp(-(log2(f02(ii,jj)./fx)*6).^2);
    end;
end;

%%%-----
function [f0,pl]=zmultiCandf0map(lfx,f0map)
%   F0 candidates extraction from time-lag representation
%   using palabolic interpolation
%   [f0,pl]=multiCand(fx,f0map)
%   lfx  : logarithmically linear frequency axis
%   f0map : time-frequency representation
%   output parameters
%   f0  : fundamental frequency (Hz)
%   pl  : peak level

%   Designed and coded by Hideki Kawahara
%   30/March/2004
%   16/June/2004 Peak picking first
%   17/June/2004 Peak selection taking into interaction account
%   26/June/2004 converted to f0map manipulation (vanilla version)

fx=log2(lfx)-log2(lfx(1));

[nr,nc]=size(f0map);
imm=diff([f0map(1,:);f0map]).*(diff([f0map;f0map(end,:)]));
dlag=diff([f0map(1,:);f0map]);
f0mapbak=f0map;
f0=zeros(nc,3);
pl=zeros(nc,3);
for ii=1:nc
    ix=find((imm(:,ii)<0)&(dlag(:,ii)>0));
    if length(ix)>0
        [mx,mxp]=max(f0map(ix,ii));
        [pl(ii,1),pos]=zParabolicInterp2(f0mapbak((-1:1)+ix(mxp),ii),ix(mxp));
    else
        pos=round(nr/4);
        pl(ii,1)=0;
    end;
    f0(ii,1)=pos*fx(2);
    if length(ix)>1
        f0map(ix(mxp),ii)=f0map(ix(mxp),ii)*0;
        [mx,mxp]=max(f0map(ix,ii));
        [pl(ii,2),pos]=zParabolicInterp2(f0mapbak((-1:1)+ix(mxp),ii),ix(mxp));
        f0(ii,2)=pos*fx(2);
        if length(ix)>2
            f0map(ix(mxp),ii)=f0map(ix(mxp),ii)*0;
            [mx,mxp]=max(f0map(ix,ii));
            [pl(ii,3),pos]=zParabolicInterp2(f0mapbak((-1:1)+ix(mxp),ii),ix(mxp));
            f0(ii,3)=pos*fx(2);
        else
            f0(ii,3)=f0(ii,2);pl(ii,3)=pl(ii,2);
        end;
    else
        f0(ii,2)=f0(ii,1);pl(ii,2)=pl(ii,1);
    end;
end;
f0=lfx(1)*2.0.^(f0);

function [val,pos]=zParabolicInterp2(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 pws=zVpowercalc(x,fs,wtc,shiftm,fc)
%   pws=Vpowercalc(x,fs,wtc,shiftm,fc)
%   x   : waveform
%   fs  : sampling frequency (Hz)
%   wtc : window time constatnt (ms)
%   shifrm  : frame update interval (ms)
%   fc  : LPF cut-off frequency (Hz)

%---- window design for pwer smoothing
t=(0:1/fs:wtc*5/1000);
w=exp(-t/(wtc/1000));
w=w-w(end);
w=w/sum(w);

%----- window for preprocesing LPF
lw=round(fs/fc*2);
b=fir1(lw-1,2*fc/fs);
nn=length(x);
x=fftfilt(b,[x(:);zeros(lw,1)]);
x=x((1:nn)+round(lw/2)-1);

yf=fftfilt(w,x.^2);
yb=fftfilt(w,x(end:-1:1).^2);
yb=yb(end:-1:1);
y=min(yf,yb);
nn=length(x);

pws=interp1((0:nn-1)/fs*1000,y,0:shiftm:(nn-1)/fs*1000);

%%%-----
function [f0r,ecr,ac1]=zrefineF06m(x,fs,f0raw,fftl,eta,nhmx,shiftm,nl,nu)
%	F0 estimation refinement
%	[f0r,ecr]=refineF06m(x,fs,f0raw,fftl,nhmx,shiftm,nl,nu)
%		x		: input waveform
%		fs		: sampling frequency (Hz)
%		f0raw	: F0 candidate (Hz)
%		fftl	: FFT length
%		eta		: temporal stretch factor
%		nhmx	: highest harmonic number
%		shiftm	: frame shift period (ms)
%		nl		: lower frame number
%		nu		: uppter frame number
%
%	Example of usage (with STRAIGHT)
%
%	global xold fs f0shiftm f0raw
%
%	dn=floor(fs/(800*3*2));
%	[f0raw,ecr]=refineF02(decimate(xold,dn),fs/dn,f0raw,512,1.1,3,f0shiftm,1,length(f0raw));

%	Designed and coded by Hideki Kawahara
%	28/July/1999
%	29/July/1999 test version using power weighting
%	30/July/1999 GcBs is added (bug fix)
%	07/August/1999 small bug fix
%   07/Dec./2002 wqitbar was added
%   30/June/2004 debug

f0raw=f0raw(:)';
f0i=f0raw;
f0i(f0i==0)=f0i(f0i==0)+160;
fx=(0:fftl-1)/fftl;
fax=(0:fftl-1)/fftl*fs;
xlinms=length(x)/fs*1000;
%nfr=min(length(f0i),round(xlinms/shiftm));
nfr=length(f0i); % 07/August/1999 

shiftl=shiftm/1000*fs;
fmx=zeros(nfr,nhmx);
vx=zeros(nfr,nhmx);
fvi=zeros(nu-nl+1,fftl/2+1);
vfs=zeros(nfr,nhmx);
avi=fvi*0;
lenx=length(x);
x=[zeros(fftl,1); x(:) ; zeros(fftl,1)]';

ec1=cos(2*pi*(0:fftl-1)/fftl); % first auto correlation basis function
ac1=f0raw*0;

tt=((1:fftl)-fftl/2)/fs;
th=(0:fftl-1)/fftl*2*pi;
rr=exp(-i*th);

f0t=100;
w1=max(0,1-abs(tt'*f0t/eta));
w1=w1(w1>0);
wg=exp(-pi*(tt*f0t/eta).^2);
wgg=(wg(abs(wg)>0.0002));
wo=fftfilt(wgg,[w1; zeros(length(wgg),1)])';

xo=(0:length(wo)-1)/(length(wo)-1);
nlo=length(wo)-1;
xi=0:1/nlo*200/100:1;
wa=interp1(xo,wo,xi,'*linear');

if nl*nu <0
	nl=1;
	nu=nfr;
end;

bx=1:fftl/2+1;
pif=zeros(fftl/2+1,nfr);
dpif=zeros(fftl/2+1,nfr);
pwm=zeros(fftl/2+1,nfr);

%%%hpg=waitbar(0,'F0 refinement using F0 adaptive analysis'); % 07/Dec./2002 by H.K.
for kk=nl:nu
		if f0i(kk) < 40
		f0i(kk)=40;
	end;
	nf0=fftl/(fs/f0i(kk));
	f0t=f0i(kk);
	xi=0:1/nlo*f0t/100:1;
	wa=interp1(xo,wo,xi,'*linear');
	wal=length(wa);
	bb=1:wal;
	bias=round(fftl-wal/2+(kk-1)*shiftl);
	dcl=mean(x(bb+bias));
    txm1=x(bb+bias-1);
    tx0=x(bb+bias);
    txp1=x(bb+bias+1);
    if (sum(abs(txm1))<1e-20)|(sum(abs(txm1))<1e-20)|(sum(abs(txm1))<1e-20)|(sum(abs(txm1))<1e-20)
        xtmp=x+randn(size(x)); % this if clause is a bug fix. 11/Jan./2005
        dcl=mean(xtmp(bb+bias));
        txm1=xtmp(bb+bias-1);
        tx0=xtmp(bb+bias);
        txp1=xtmp(bb+bias+1);
    end;
    %if sum(abs(txm1))==0;txm1=randn(size(txm1));end;
    %if sum(abs(tx0))==0;tx0=randn(size(tx0));end;
    %if sum(abs(txp1))==0;txp1=randn(size(txp1));end;
	%ff0=fft((x(bb+bias-1)-dcl).*wa,fftl);
	%ff1=fft((x(bb+bias)-dcl).*wa,fftl);
	%ff2=fft((x(bb+bias+1)-dcl).*wa,fftl);
	ff0=fft((txm1-dcl).*wa,fftl);
	ff1=fft((tx0-dcl).*wa,fftl);
	ff2=fft((txp1-dcl).*wa,fftl);
    ff0(ff0==0)=ff0(ff0==0)+0.000000001;
    ff1(ff1==0)=ff1(ff1==0)+0.000000001;
    ff2(ff2==0)=ff2(ff2==0)+0.000000001;
	fd=ff2.*rr-ff1;
	fd0=ff1.*rr-ff0;
	crf=fax+(real(ff1).*imag(fd)-imag(ff1).*real(fd))./(abs(ff1).^2)*fs/pi/2;
	crf0=fax+(real(ff0).*imag(fd0)-imag(ff0).*real(fd0))./(abs(ff0).^2)*fs/pi/2;
	pif(:,kk)=crf(bx)'*2*pi;
	dpif(:,kk)=(crf(bx)-crf0(bx))'*2*pi;
	pwm(:,kk)=abs(ff1(bx)'); % 29/July/1999
    %%%waitbar((kk-nl)/(nu-nl));% ,hpg) % 07/Dec./2002 by H.K.
    ac1(kk)=sum(abs(ff1).^2.0.*ec1)/sum(abs(ff1).^2);
%	keyboard;
end;
%%%close(hpg);
slp=([pif(2:fftl/2+1,:);pif(fftl/2+1,:)]-pif)/(fs/fftl*2*pi);
dslp=([dpif(2:fftl/2+1,:);dpif(fftl/2+1,:)]-dpif)/(fs/fftl*2*pi)*fs;
mmp=slp*0;

%[c1,c2]=znrmlcf2(1);
[c1,c2]=znrmlcf3(shiftm);
fxx=((0:fftl/2)+0.5)/fftl*fs*2*pi;

%--- calculation of relative noise level

%%%hpg=waitbar(0,'P/N calculation'); % 07/Dec./2002 by H.K.
for ii=1:fftl/2+1;
	c2=c2*(fxx(ii)/2/pi)^2;
	mmp(ii,:)=(dslp(ii,:)/sqrt(c2)).^2+(slp(ii,:)/sqrt(c1)).^2;
    %%%if rem(ii,10)==0;waitbar(ii/(fftl/2+1));end;  % 07/Dec./2002 by H.K.
end;
%%%close(hpg); % 07/Dec./2002 by H.K.

%--- Temporal smoothing

%sml=round(6*fs/1000/2)*2+1; % 12 ms, and odd number
%sml=round(4*fs/1000/2)*2+1; % 8 ms, and odd number
sml=round(1.5*fs/1000/2/shiftm)*2+1; % 3 ms, and odd number
smb=round((sml-1)/2); % bias due to filtering

%%%hpg=waitbar(0,'P/N smoothing'); % 07/Dec./2002 by H.K.
%This smoothing is modified (30 Nov. 2000).