-
Notifications
You must be signed in to change notification settings - Fork 29
/
refineF06.m
203 lines (175 loc) · 6.36 KB
/
refineF06.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
function [f0r,ecr]=refineF06(x,fs,f0raw,fftl,eta,nhmx,shiftm,nl,nu,imgi);
% F0 estimation refinement
% [f0r,ecr]=refineF06(x,fs,f0raw,fftl,nhmx,shiftm,nl,nu,imgi)
% 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
% imgi : display indicator, 1: display on (default), 0: off
%
% 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
% 13.May/2005 minor vulnerability fix
% 10/Aug./2005 modified by Takahashi on waitbar
% 10/Sept./2005 modified by Kawahara on waitbar
% 16/Sept./2005 minor bug fix
% 26/Sept./2005 bug fix
if nargin==9; imgi=1; end;
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)]';
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);
rmsValue = std(x); % 26/Sept./2005 by HK
if imgi==1; hpg=waitbar(0,'F0 refinement using F0 adaptive analysis'); end; % 07/Dec./2002 by H.K.%10/Aug./2005
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);
if std(x(bb+bias))*std(x(bb+bias-1))*std(x(bb+bias+1)) == 0 % 26/Sept./2005 by HK
x(bb+bias) = randn(length(bb),1)*rmsValue/100000;
end;
dcl=mean(x(bb+bias));
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);
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
if imgi==1; waitbar((kk-nl)/(nu-nl)); end; % ,hpg) % 07/Dec./2002 by H.K.%10/Aug./2005
% keyboard;
end;
if imgi==1; close(hpg); end;
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]=znrmlcf2(shiftm);
fxx=((0:fftl/2)+0.5)/fftl*fs*2*pi;
%--- calculation of relative noise level
if imgi==1; hpg=waitbar(0,'P/N calculation'); end; % 07/Dec./2002 by H.K.%10/Aug./2005
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 imgi==1 & rem(ii,10)==0;waitbar(ii/(fftl/2+1));end; % 07/Dec./2002 by H.K.%10/Aug./2005
end;
if imgi==1; close(hpg); end; % 07/Dec./2002 by H.K.%10/Aug./2005
%--- 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
if imgi==1; hpg=waitbar(0,'P/N smoothing'); end; % 07/Dec./2002 by H.K.%10/Aug./2005
%This smoothing is modified (30 Nov. 2000).
smmp=fftfilt((hanning(sml).^2)/sum((hanning(sml).^2)), ...
[mmp zeros(fftl/2+1,sml*2)]'+max(max(mmp((~isnan(mmp))&(mmp<Inf))))*0.0000001)'; % 26/Sept./2005 by H.K.
if imgi==1; waitbar(0.2); end;%10/Aug./2005
smmp=1.0./fftfilt(hanning(sml)/sum(hanning(sml)),1.0./smmp')';
if imgi==1; waitbar(0.4); end;
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)]'+0.00001)';
if imgi==1; waitbar(0.6); end;%10/Aug./2005
spfm=fftfilt(hanning(sml)/sum(hanning(sml)),[pwm.*pif zeros(fftl/2+1,sml*2)]'+0.00001)';
if imgi==1; waitbar(0.8); end;%10/Aug./2005
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;
if imgi==1; waitbar(1); end;%10/Aug./2005
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);
if imgi==1; close(hpg); end;%10/Aug./2005
%keyboard;
%--------------------
function [c1,c2]=znrmlcf2(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;