multicuef0.m

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

%zz=0:1/n:3;
%hh=[diff(zGcBs(zz,0)) 0]*n;
%c1=sum((zz.*hh).^2)/n;
%c2=sum((2*pi*zz.^2.*hh).^2)/n;

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

function p=zGcBs(x,k)

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


%--------------------------------------------
function smap=zsmoothmapB(map,fs,f0floor,nvo,mu,mlim,pex)

[nvc,mm]=size(map);
%mu=0.4;
t0=1/f0floor;
lmx=round(6*t0*fs*mu);
wl=2^ceil(log(lmx)/log(2));
gent=((1:wl)-wl/2)/fs;

smap=map;
mpv=1;
zt=0*gent;
iiv=1:mm;
for ii=1:nvc
    t=gent*mpv; %t0*mu/mpv*1000
    t=t(abs(t)<3.5*mu*t0);
    wbias=round((length(t)-1)/2);
    wd1=exp(-pi*(t/(t0*(1-pex))/mu).^2);
    wd2=exp(-pi*(t/(t0*(1+pex))/mu).^2);
    wd1=wd1/sum(wd1);
    wd2=wd2/sum(wd2);
    tm=fftfilt(wd1,[map(ii,:) zt]);
    tm=fftfilt(wd2,[1.0./tm(iiv+wbias) zt]);
    smap(ii,:)=1.0./tm(iiv+wbias);
    if t0*mu/mpv*1000 > mlim
        mpv=mpv*(2.0^(1/nvo));
    end;
end;

%--------------------------------------------
function [ff,vv,df]=zfixpfreq2(fxx,pif2,mmp,dfv)

nn=length(fxx);
iix=(1:nn)';
cd1=pif2-fxx;
cd2=[diff(cd1);cd1(nn)-cd1(nn-1)];
cdd1=[cd1(2:nn);cd1(nn)];
fp=(cd1.*cdd1<0).*(cd2<0);
ixx=iix(fp>0);
ff=pif2(ixx)+(pif2(ixx+1)-pif2(ixx)).*cd1(ixx)./(cd1(ixx)-cdd1(ixx));
%vv=mmp(ixx);
vv=mmp(ixx)+(mmp(ixx+1)-mmp(ixx)).*(ff-fxx(ixx))./(fxx(ixx+1)-fxx(ixx));
df=dfv(ixx)+(dfv(ixx+1)-dfv(ixx)).*(ff-fxx(ixx))./(fxx(ixx+1)-fxx(ixx));

%--------------------------------------------
function [ff,vv,df,aa]=zfixpfreq3(fxx,pif2,mmp,dfv,pm)

aav=abs(pm);
nn=length(fxx);
iix=(1:nn)';
cd1=pif2-fxx;
cd2=[diff(cd1);cd1(nn)-cd1(nn-1)];
cdd1=[cd1(2:nn);cd1(nn)];
fp=(cd1.*cdd1<0).*(cd2<0);
ixx=iix(fp>0);
ff=pif2(ixx)+(pif2(ixx+1)-pif2(ixx)).*cd1(ixx)./(cd1(ixx)-cdd1(ixx));
%vv=mmp(ixx);
vv=mmp(ixx)+(mmp(ixx+1)-mmp(ixx)).*(ff-fxx(ixx))./(fxx(ixx+1)-fxx(ixx));
df=dfv(ixx)+(dfv(ixx+1)-dfv(ixx)).*(ff-fxx(ixx))./(fxx(ixx+1)-fxx(ixx));
aa=aav(ixx)+(aav(ixx+1)-aav(ixx)).*(ff-fxx(ixx))./(fxx(ixx+1)-fxx(ixx));

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

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

%--------------------------------------------
function x=cleaninglownoise(x,fs,f0floor);

flm=50;
flp=round(fs*flm/1000);
nn=length(x);
wlp=fir1(flp*2,f0floor/(fs/2));
wlp(flp+1)=wlp(flp+1)-1;
wlp=-wlp;

tx=[x(:)' zeros(1,2*length(wlp))];
ttx=fftfilt(wlp,tx);
x=ttx((1:nn)+flp);

return;

%%%---------
function pm=zmultanalytFineCSPB(x,fs,f0floor,nvc,nvo,mu,mlt);

%       Dual waveleta analysis using cardinal spline manipulation
%               pm=multanalytFineCSPB(x,fs,f0floor,nvc,nvo,mu,mlt)
%       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
%				mlt		: harmonic ID#
%       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
%       07/Dec./2002 waitbar was added

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)];
txx=tx;
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;
%fs
%mu=1.0;
%%%hpg=waitbar(0,['wavelet analysis for initial F0 and P/N estimation with HM#:' num2str(mlt)]); % 07/Dec./2002 by H.K.
for ii=1:nvc
  tb=gent*mpv;
  t=tb(abs(tb)<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.00001);
  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);
%  ampp=abs(pmtmp1(wbias+1:wbias+nx));
%  txx=[x./ampp,zeros(1,wl)];
%  txx=tx;
%  pmtmp1=fftfilt(wwd,txx);  
%  disp(['ii= ' num2str(ii) '  sum=' num2str(sum(abs(wwd)),6)]);
  pm(ii,:)=pmtmp1(wbias+1:wbias+nx)*sqrt(mpv);
  mpv=mpv*(2.0^(1/nvo));
  %%%waitbar(ii/nvc);%,hpg); % 07/Dec./2002 by H.K.
%  keyboard;
end;
%%%close(hpg); % 07/Dec./2002 by H.K.
%[nn,mm]=size(pm);
%pm=pm(:,1:mlt:mm);

%%%----
function x=zdB(y)
x=20*log10(y);

%%%-----
function [val,pos]=zmultiCandIF(f0v,vrv)
%   [val,pos]=multiCandIF(f0v,vrv)
%   F0 candidates based on instantaneous frequency
%   fixed points
%   f0v : fixed point frequencies (Hz)
%   vrv : fixed point N/C (ratio)

%   by Hideki Kawahara
%   23/June/2004

[nr,nc]=size(f0v);
[nr2,nc2]=size(vrv);
if (nr~=nr2) | (nc~=nc2);val=[];pos=[];return;end;
vrvdb=-zdBpower(vrv);

mxfq=100000;
val=zeros(nc,3);
pos=ones(nc,3);
for ii=1:nc
    f=f0v(:,ii)';
    v=vrvdb(:,ii)';
    v=v(f<mxfq);
    f=f(f<mxfq);
    [mx,mxp]=max(v);
    if length(mxp)>0
        pos(ii,1)=f(mxp);
        val(ii,1)=v(mxp);
    else
        pos(ii,1)=1;
        val(ii,1)=1;
    end;
    if length(f)>1
        v(mxp)=v(mxp)*0-50;
        [mx,mxp]=max(v);
        pos(ii,2)=f(mxp);
        val(ii,2)=v(mxp);
        if length(f)>2
            v(mxp)=v(mxp)*0-50;
            [mx,mxp]=max(v);
            pos(ii,3)=f(mxp);
            val(ii,3)=v(mxp);
        else
            pos(ii,3)=pos(ii,2);val(ii,3)=val(ii,2);
        end;
    else
        pos(ii,2)=pos(ii,1);val(ii,2)=val(ii,1);
    end;
end;

%%%------
function y=zdBpower(x)
y=10*log10(x);

%%%------
function [y,ind,fq]=zremoveACinduction(x,fs,pos);
%   [y,ind,fq]=removeACinduction(x,fs,pos);
%   Function to remove AC induction
%   x   : input speech signal 
%   fs  : sampling frequency (Hz)
%   pos : Locations of Top-three F0 candidates (Hz)
%   Output parameter
%   y   : speech signal without AC induction
%   ind : 1 indicates AC induction was detected
%   fq  : frequency of AC induction

%   Designed and coded by Hideki Kawahawra, 
%   24/June/2004

x=x(:);
ind=0;
f=pos(:);
h50=sum(abs(f-50)<5)/sum(f>0);
h60=sum(abs(f-60)<5)/sum(f>0);
if (h50<0.2) & (h60<0.2);y=x;fq=0;return;end;
ind=1;
if h50>h60
    fq=50;
else
    fq=60;
end;
tx=(1:length(x))'/fs;
fqv=([-0.3:0.025:0.3]+fq);
txv=tx*fqv;
fk=x'*exp(-j*2*pi*txv)/length(x);
[fkm,ix]=max(abs(fk));
fq=fqv(ix);
y=x-2*real(fk(ix)*exp(j*2*pi*fq*tx));

%%%----
function [lagspec,lx]=zlagspectestnormal(x,fs,stp,edp,shiftm,wtlm,ndiv,wflf,pc,amp,h)
%   Lag spectrogram for F0 extraction
%   [lagspec,lx]=lagspectestnormal(x,fs,stp,edp,shiftm,wtlm,ndiv,wflf,pc,amp,h)
%   x   : waveform
%   fs  : sampling frequency (Hz)
%   stp : starting position (ms)
%   edp : end position (ms)
%   shiftm  : frame shift for analysis (ms)
%   wtlm    : time window length (ms)
%   ndiv    : number of segment in the frequency domain
%   wflf    : frequency domain window length (Hz)
%   pc      : power exponent for nonlinearity
%   amp : amount of lag window compensation

%   h       : handle for graph

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

nftm=floor((edp-stp)/shiftm);
pm=stp;
[acc,abase,fx,lx]=ztestspecspecnormal(x,fs,pm,wtlm,ndiv,wflf,pc,amp);
nlx=length(lx);
lagspec=zeros(nlx,nftm);
for ii=1:nftm
    pmmul=stp+(ii-1)*shiftm;
    [acc,abase,fx,lx]=ztestspecspecnormal(x,fs,pmmul,wtlm,ndiv,wflf,pc,amp);%keyboard;
    lagspec(:,ii)=mean(acc')'/mean(acc(1,:));
end;
if h>0
    figure(h);
    imagesc([stp edp],[0 max(lx)]*1000,max(0,lagspec));
    axis('xy')
    axis([stp edp 0 40]);
    set(gca,'fontsize',16);
    xlabel('time (ms)')
    ylabel('lag (ms)')
    title(['wtl=' num2str(wtlm) 'ms ndiv=' num2str(ndiv) ' wfl=' num2str(wflf) 'Hz PC=' num2str(pc) ...
            ' fs=' num2str(fs) 'Hz amp=' num2str(amp)]);
    drawnow;
end;

%%%------
function [acc,abase,fx,lx]=ztestspecspecsimple(x,fs,pm,wtlm,ndiv,wflf,pc,amp);
%   Modified auto correlation 
%   [acc,abase,fx,lx]=testspecspecsimple(x,fs,pm,wtlm,ndiv,wflf,pc,amp);
%   input parameters
%   x   : signal to be analyzed
%   fs  : sampling frequency (Hz)
%   pm   : position to be tested (ms)
%   wtlm    : time window length (ms)
%   ndiv    : number of division on frequency axis
%   wflf    : frequency window length (hz)
%   pc  : power exponent
%   amp : amount of lag window compensation
%   output parameters
%   acc : spectrogram on frequency axis
%       : (periodicity gram on local frequency area)
%   fx  : frequency axis
%   lx  : lag axis

%  Test program for spectrum check
%   by Hideki Kawahara 27 March 2004
%   29/March/2004 streamlined
%   12/June/2004 Bias term removed
%   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
abase=abs(fft(xt.*blackman(wtlmso),fftl));
mxab=max(abase);
ac=ifft(abase.^2);
npw=real(fft(ac.*lagw'));
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 [acc,abase,fx,lx]=ztestspecspecnormal(x,fs,pm,wtlm,ndiv,wflf,pc,amp);
%   Modified auto correlation 
%   [acc,abase,fx,lx]=testspecspecnormal(x,fs,pm,wtlm,ndiv,wflf,pc,amp);
%   input parameters
%   x   : signal to be analyzed
%   fs  : sampling frequency (Hz)
%   pm   : position to be tested (ms)
%   wtlm    : time window length (ms)
%   ndiv    : number of division on frequency axis
%   wflf    : frequency window length (hz)
%   pc  : power exponent
%   amp : amount of lag window compensation
%   output parameters
%   acc : spectrogram on frequency axis
%       : (periodicity gram on local frequency area)
%   fx  : frequency axis
%   lx  : lag axis

%  Test program for spectrum check
%   by Hideki Kawahara 27 March 2004
%   29/March/2004 streamlined
%   12/June/2004 Bias term removed