multicuef0.m

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

%thsilence=3; % for silence
%thvuv=0.6; % for V/UV
%f00=f0raw2;
%f00(pwsdb(1:length(f00))<(noiselevel+thsilence))=f00(pwsdb(1:length(f00))<(noiselevel+thsilence))*0;
%f0raw3=f00;
%f0raw3((zdB(ecr)/100+ac1)<thvuv)=f0raw3((zdB(ecr)/100+ac1)<thvuv)*0;
%----- new V/UV decision routine 15/Aug./2004
auxouts.BackgroundNoiselevel=noiselevel;
vuv=zvuvdecision4(f0raw2,auxouts);
nnll=min(length(f0raw2),length(vuv));
f0raw3=f0raw2(1:nnll).*vuv(1:nnll);
if imgi==1
    semilogy(f0raw3,'k');hold off
    title('F0 estimates,  cyan:initial, greeen:fine-tuned, black:voiced part')
    xlabel('time (ms)');ylabel('frequency (Hz)');
end;
%vuv=zdB(ecr)/100+ac1;
f0raw=f0raw2(1:nnll);vuv=vuv(1:nnll);
auxouts.F0candidatesByIF=pos;
auxouts.CNofcandidatesByIF=val;
auxouts.F0candidatesByAC=f02;
auxouts.ACofcandidatesByAC=pl2;
auxouts.F0candidatesByMix=f0cand;
auxouts.RELofcandidatesByMix=relv;
auxouts.RefinedCN=ecr;
auxouts.FirstAutoCorrelation=ac1;
auxouts.F0initialEstimate=f0raw0;
auxouts.BackgroundNoiselevel=noiselevel;
auxouts.InstantaneousPower=pwsdb;
auxouts.RefinedF0estimates=f0raw;
auxouts.VUVindicator=vuv;
if imgi==1; displaysummary(auxouts,f0floor,f0ceil); end;
switch nargout
    case 0
        f0raw=auxouts;eval(['help ' mfilename]);
    case 1
        f0raw=f0raw3;
    case 2
    case {3,4}
    otherwise
        eval(['help ' mfilename]);
        return;
end

%------
function prm=zsetdefaultparams;
prm.F0searchLowerBound=40; % f0floor
prm.F0searchUpperBound=800; % f0ceil
prm.F0frameUpdateInterval=1; % shiftm % F0 calculation interval (ms)
prm.NofChannelsInOctave=24; % nvo=24; % Number of channels in one octave
prm.IFWindowStretch=1.2; % mu=1.2; % window stretch from isometric window
prm.DisplayPlots=0; % imgi=1; % image display indicator (1: display image)
prm.IFsmoothingLengthRelToFc=1; %  smp=1; % smoothing length relative to fc (ratio)
prm.IFminimumSmoothingLength=5; %  minm=5; % minimum smoothing length (ms)
prm.IFexponentForNonlinearSum=0.5; % pc=0.5; % exponent to represent nonlinear summation
prm.IFnumberOfHarmonicForInitialEstimate=1; % nc=1; % number of harmonic component to use (1,2,3)
prm.TimeConstantForPowerCalculation=10;  % tcpower=10; % time constant for power calculation (ms)
prm.ACtimeWindowLength=60; % Time window length for Autocorrelation based method (ms)
prm.ACnumberOfFrequencySegments=8; % for Autocorrelation method
prm.ACfrequencyDomainWindowWidth=2200; % for Autocorrelation method (Hz)
prm.ACpowerExponentForNonlinearity=0.5; % for Autocorrelation method
prm.ACamplitudeCompensationInShortLag=1.6; %2.2; % for Autocorrelation method (ratio) 23/July/2004
prm.ACexponentForACdistance=4; % Nonlinear distance measure for post processing
prm.AClagSmoothingLength=0.0001; % Lag smoothing length for post processing (s) !! 0.01 to 0.0001 23/July/2004
prm.ACtemporalSmoothingLength=20; % Temporal smoothing length for post processing (ms)
prm.ThresholdForSilence=3; % for silence decision above average noise level (dB)
prm.ThresholdForVUV=0.6; % for V/UV decision based on first autocorrelation and C/N
prm.WeightForAutocorrelationMap=1; % weight for combining maps (Autocorrelation)
prm.WeightForInstantaneousFqMap=1; % weight for combining maps (Instantaneous Frequency)
prm.VUVthresholdOfAC1=-0.1; % First autocorrelation thershould for VUV in segment search
prm.SDforNormalizeMixingDistance=0.3; % Normalization factor for mixing F0 distance (octave)
prm.SDforTrackingNormalization=0.2;
prm.MaxumumPermissibleOctaveJump=0.4; 
prm.ThresholdToStartSearch=0.3;
prm.ThresholdToQuitSearch=0.35;
prm.ThresholdForReliableRegion=0.25;
prm.WhoAmI=mfilename;
return;

%%%------------
function oki=displaysummary(f,f0floor,f0ceil)
oki=1;
nn=length(f.RefinedF0estimates);
figure
subplot(211);
semilogy(f.F0initialEstimate,'c');grid on;
hold on;
semilogy(f.RefinedF0estimates.*f.VUVindicator,'b');grid on;
set(gca,'fontsize',14);
xlabel('time (ms)');
ylabel('frequency (Hz)');
axis([1 nn f0floor f0ceil]);
subplot(212);
plot(f.RELofcandidatesByMix,'.');grid on;
set(gca,'fontsize',14);
xlabel('time (ms)');
ylabel('relative periodicity');
axis([1 nn 0 1]);
drawnow

%%%------------
function [f0v,vrv,dfv,nf,aav]=zfixpF0VexMltpBG4(x,fs,f0floor,nvc,nvo,mu,imgi,shiftm,smp,minm,pc,nc)

%	Fixed point analysis to extract F0
%	[f0v,vrv,dfv,nf]=fixpF0VexMltpBG4(x,fs,f0floor,nvc,nvo,mu,imgi,shiftm,smp,minm,pc,nc)
%	x	: input signal
%	fs	: sampling frequency (Hz)
%	f0floor	: lowest frequency for F0 search
%	nvc	: total number of filter channels
%	nvo	: number of channels per octave
%	mu	: temporal stretching factor
%	imgi	: image display indicator (1: display image)
%	shiftm	: frame shift in ms
%	smp	: smoothing length relative to fc (ratio)
%	minm	: minimum smoothing length (ms)
%	pc	: exponent to represent nonlinear summation
%	nc	: number of harmonic component to use (1,2,3)

%	Designed and coded by Hideki Kawahara
%	28/March/1999
%	04/April/1999 revised to multi component version
%	07/April/1999 bi-reciprocal smoothing for multi component compensation
%	01/May/1999 first derivative of Amplitude is taken into account
%	17/Dec./2000 display bug fix
%	19/Sep./2002 bug fix (mu information was discarded.)
%   07/Dec./2002 waitbar was added

%f0floor=40;
%nvo=12;
%nvc=52;
%mu=1.1;

x=cleaninglownoise(x,fs,f0floor);

fxx=f0floor*2.0.^((0:nvc-1)/nvo)';
fxh=max(fxx);

dn=max(1,floor(fs/(fxh*6.3)));

if nc>2
	pm3=zmultanalytFineCSPB(decimate(x,dn),fs/dn,f0floor,nvc,nvo,mu,3); % error crrect 2002.9.19 (mu was fixed 1.1)
	pif3=zwvlt2ifq(pm3,fs/dn);
	[nn,mm]=size(pif3);
	pif3=pif3(:,1:3:mm);
	pm3=pm3(:,1:3:mm);
end;

if nc>1
	pm2=zmultanalytFineCSPB(decimate(x,dn),fs/dn,f0floor,nvc,nvo,mu,2);% error crrect 2002.9.19(mu was fixed 1.1)
	pif2=zwvlt2ifq(pm2,fs/dn);
	[nn,mm]=size(pif2);
	pif2=pif2(:,1:3:mm);
	pm2=pm2(:,1:3:mm);
end;

pm1=zmultanalytFineCSPB(decimate(x,dn*3),fs/(dn*3),f0floor,nvc,nvo,mu,1);% error crrect 2002.9.19(mu was fixed 1.1)
%%%% safe guard added on 15/Jan./2003
mxpm1=max(max(abs(pm1)));
eeps=mxpm1/10000000;
pm1(pm1==0)=pm1(pm1==0)+eeps;
%%%% safe guard end
pif1=zwvlt2ifq(pm1,fs/(dn*3));
%keyboard;

[nn,mm1]=size(pif1);
mm=mm1;
if nc>1
	[nn,mm2]=size(pif2);
	mm=min(mm1,mm2);
end;

if nc>2
	[nn,mm3]=size(pif3);
	mm=min([mm1 mm2 mm3]);
end;

if nc == 2
	for ii=1:mm
		pif2(:,ii)=(pif1(:,ii).*(abs(pm1(:,ii))).^pc ...
			+pif2(:,ii)/2.*(abs(pm2(:,ii))).^pc )...
			./((abs(pm1(:,ii))).^pc+(abs(pm2(:,ii))).^pc);
	end;
end;
if nc == 3 
	for ii=1:mm
		pif2(:,ii)=(pif1(:,ii).*(abs(pm1(:,ii))).^pc ...
			+pif2(:,ii)/2.*(abs(pm2(:,ii))).^pc ...
			+pif3(:,ii)/3.*(abs(pm3(:,ii))).^pc )... 
			./((abs(pm1(:,ii))).^pc+(abs(pm2(:,ii))).^pc+(abs(pm3(:,ii))).^pc);
	end;
end;
if nc == 1
	pif2=pif1;
end;

		  
%pif2=zwvlt2ifq(pm,fs/dn)*2*pi;
pif2=pif2*2*pi;
dn=dn*3;

[slp,pbl]=zifq2gpm2(pif2,f0floor,nvo);
[nn,mm]=size(pif2);
dpif=[pif2(:,2:mm)-pif2(:,1:mm-1)]*fs/dn;
dpif(:,mm)=dpif(:,mm-1);
[dslp,dpbl]=zifq2gpm2(dpif,f0floor,nvo);

damp=[abs(pm1(:,2:mm))-abs(pm1(:,1:mm-1))]*fs/dn;
damp(:,mm)=damp(:,mm-1);
damp=damp./abs(pm1);

%[c1,c2]=znormwght(1000);
fxx=f0floor*2.0.^((0:nn-1)/nvo)'*2*pi;
mmp=0*dslp;
[c1,c2b]=znrmlcf2(1);
%%%hpg=waitbar(0,'P/N map calculation'); % 07/Dec./2002 by H.K.
for ii=1:nn
%	[c1,c2]=znrmlcf2(fxx(ii)/2/pi); % This is OK, but the next Eq is much faster.
	c2=c2b*(fxx(ii)/2/pi)^2;
	cff=damp(ii,:)/fxx(ii)*2*pi*0;
	mmp(ii,:)=(dslp(ii,:)./(1+cff.^2)/sqrt(c2)).^2+(slp(ii,:)./sqrt(1+cff.^2)/sqrt(c1)).^2;
    %%%waitbar(ii/nn);%,hpg); % 07/Dec./2002 by H.K.
end;
%%%close(hpg);

if smp~=0
	smap=zsmoothmapB(mmp,fs/dn,f0floor,nvo,smp,minm,0.4);
else
	smap=mmp;
end;

fixpp=zeros(round(nn/3),mm);
fixvv=fixpp+100000000;
fixdf=fixpp+100000000;
fixav=fixpp+1000000000;
nf=zeros(1,mm);
%%%hpg=waitbar(0,'Fixed pints calculation'); % 07/Dec./2002 by H.K.
for ii=1:mm
	[ff,vv,df,aa]=zfixpfreq3(fxx,pif2(:,ii),smap(:,ii),dpif(:,ii)/2/pi,pm1(:,ii));
	kk=length(ff);
	fixpp(1:kk,ii)=ff;
	fixvv(1:kk,ii)=vv;
	fixdf(1:kk,ii)=df;
	fixav(1:kk,ii)=aa;
	nf(ii)=kk;
    %%%if rem(ii,10)==0; waitbar(ii/mm); end;% 07/Dec./2002 by H.K.
end;
%%%close(hpg); % 07/Dec./2002 by H.K.
fixpp(fixpp==0)=fixpp(fixpp==0)+1000000;

%keyboard
%[vvm,ivv]=min(fixvv);
%
%for ii=1:mm
%	ff00(ii)=fixpp(ivv(ii),ii);
%	esgm(ii)=fixvv(ivv(ii),ii);
%end;
np=max(nf);
f0v=fixpp(1:np,round(1:shiftm/dn*fs/1000:mm))/2/pi;
vrv=fixvv(1:np,round(1:shiftm/dn*fs/1000:mm));
dfv=fixdf(1:np,round(1:shiftm/dn*fs/1000:mm));
aav=fixav(1:np,round(1:shiftm/dn*fs/1000:mm));
nf=nf(round(1:shiftm/dn*fs/1000:mm));

if imgi==1
	okid=cnmap(fixpp,smap,fs,dn,nvo,f0floor,shiftm);
end;
%ff00=ff00(round(1:shiftm/dn*fs/1000:mm));
%esgm=sqrt(esgm(round(1:shiftm/dn*fs/1000:mm)));
%keyboard;

return;
%------------------------------------------------------------------
function okid=cnmap(fixpp,smap,fs,dn,nvo,f0floor,shiftm)

% 	This function had a bug in map axis.
%	17/Dec./2000 bug fix by Hideki Kawahara.

dt=dn/fs;
[nn,mm]=size(smap);
aa=figure;
%set(aa,'PaperPosition',[0.3 0.25 8 10.9]);
%set(aa,'Position',[30 130 520 680]);
%subplot(211);
imagesc([0 (mm-1)*dt*1000],[1 nn],zdB(smap(:,round(1:shiftm/dn*fs/1000:mm))));axis('xy')
set(gca,'fontsize',16);
hold on;
tx=((1:shiftm/dn*fs/1000:mm)-1)*dt*1000;
plot(tx,(nvo*log(fixpp(:,round(1:shiftm/dn*fs/1000:mm))/f0floor/2/pi)/log(2)+0.5)','ko');
plot(tx,(nvo*log(fixpp(:,round(1:shiftm/dn*fs/1000:mm))/f0floor/2/pi)/log(2)+0.5)','w.');
hold off
xlabel('time (ms)');
ylabel('channel #');
title('Instantaneous frequency-based fixed points')
colormap(jet);

okid=1;
return;

%------------------------------------------------------------------

function pm=zmultanalytFineCSPm(x,fs,f0floor,nvc,nvo,mu,mlt);

%       Dual waveleta analysis using cardinal spline manipulation
%               pm=multanalytFineCSP(x,fs,f0floor,nvc,nvo);
%       Input parameters 
%               
%               x       : input signal (2kHz sampling rate is sufficient.)
%               fs      : sampling frequency (Hz)
%               f0floor : lower bound for pitch search (60Hz suggested)
%               nvc     : number of total voices for wavelet analysis
%               nvo     : number of voices in an octave
%				mu		: temporal stretch factor
%       Outpur parameters
%               pm      : wavelet transform using iso-metric Gabor function
%
%       If you have any questions,  mailto:[e-mail]
%
%       Copyright (c) ATR Human Information Processing Research Labs. 1996
%       Invented and coded by Hideki Kawahara
%       30/Oct./1996

t0=1/f0floor;
lmx=round(6*t0*fs*mu);
wl=2^ceil(log(lmx)/log(2));
x=x(:)';
nx=length(x);
tx=[x,zeros(1,wl)];
gent=((1:wl)-wl/2)/fs;

%nvc=18;

wd=zeros(nvc,wl);
wd2=zeros(nvc,wl);
ym=zeros(nvc,nx);
pm=zeros(nvc,nx);
mpv=1;
%mu=1.0;
for ii=1:nvc
    t=gent*mpv;
    t=t(abs(t)<3.5*mu*t0);
    wbias=round((length(t)-1)/2);
    wd1=exp(-pi*(t/t0/mu).^2);%.*exp(i*2*pi*t/t0); 
    wd2=max(0,1-abs(t/t0/mu));
    wd2=wd2(wd2>0);
    wwd=conv(wd2,wd1);
    wwd=wwd(abs(wwd)>0.0001);
    wbias=round((length(wwd)-1)/2);
    wwd=wwd.*exp(i*2*pi*mlt*t(round((1:length(wwd))-wbias+length(t)/2))/t0);
    pmtmp1=fftfilt(wwd,tx);
    pm(ii,:)=pmtmp1(wbias+1:wbias+nx)*sqrt(mpv);
    mpv=mpv*(2.0^(1/nvo));
    %  keyboard;
end;
%[nn,mm]=size(pm);
%pm=pm(:,1:mlt:mm);

%----------------------------------------------------------------
function pif=zwvlt2ifq(pm,fs)
%	Wavelet to instantaneous frequency map
%	fqv=wvlt2ifq(pm,fs)

%	Coded by Hideki Kawahara
%	02/March/1999

[nn,mm]=size(pm);
pm=pm./(abs(pm));
pif=abs(pm(:,:)-[pm(:,1),pm(:,1:mm-1)]);
pif=fs/pi*asin(pif/2);
pif(:,1)=pif(:,2);

%----------------------------------------------------------------

function [slp,pbl]=zifq2gpm2(pif,f0floor,nvo)
%	Instantaneous frequency 2 geometric parameters
%	[slp,pbl]=ifq2gpm(pif,f0floor,nvo)
%	slp		: first order coefficient
%	pbl		: second order coefficient

%	Coded by Hideki Kawahara
%	02/March/1999

[nn,mm]=size(pif);
fx=f0floor*2.0.^([0:nn-1]/nvo)*2*pi;

c=2.0^(1/nvo);
g=[1/c/c 1/c 1;1 1 1;c*c c 1];
h=inv(g);

%slp=pif(1:nn-2,:)*h(1,1)+pif(2:nn-1,:)*h(1,2)+pif(3:nn,:)*h(1,3);
slp=((pif(2:nn-1,:)-pif(1:nn-2,:))/(1-1/c) ...
    +(pif(3:nn,:)-pif(2:nn-1,:))/(c-1))/2;
slp=[slp(1,:);slp;slp(nn-2,:)];

pbl=pif(1:nn-2,:)*h(2,1)+pif(2:nn-1,:)*h(2,2)+pif(3:nn,:)*h(2,3);
pbl=[pbl(1,:);pbl;pbl(nn-2,:)];

for ii=1:nn
	slp(ii,:)=slp(ii,:)/fx(ii);
	pbl(ii,:)=pbl(ii,:)/fx(ii);
end;

%------------------------------------------

%function [c1,c2]=znormwght(n)