0001 function [fx,tx,pv,fv]=v_fxpefac(s,fs,tinc,m,pp)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055 persistent w_u m_u v_u w_v m_v v_v dpwtdef
0056
0057 if ~numel(w_u)
0058
0059
0060
0061
0062
0063 if nargin>3 && any(m=='x')
0064 fxpefac_g;
0065 fxpefac_w;
0066 else
0067 w_u=[0.1461799 0.3269458 0.2632178 0.02331986 0.06360947 0.1767271 ]';
0068
0069 m_u=[13.38533 0.4199435 ;
0070 12.23505 0.1496836 ;
0071 12.76646 0.2581733 ;
0072 13.69822 0.6893078 ;
0073 9.804372 0.02786567 ;
0074 11.03848 0.07711229 ];
0075
0076 v_u=reshape([0.4575519 0.002619074 0.002619074 0.01262138 ;
0077 0.7547719 0.008568089 0.008568089 0.001933864 ;
0078 0.5770533 0.003561592 0.003561592 0.00527957 ;
0079 0.3576287 0.01388739 0.01388739 0.04742106 ;
0080 0.9049906 0.01033191 0.01033191 0.0001887114 ;
0081 0.637969 0.009936445 0.009936445 0.0007082946 ]',[2 2 6]);
0082
0083 w_v=[0.1391365 0.221577 0.2214025 0.1375109 0.1995124 0.08086066 ]';
0084
0085 m_v=[15.36667 0.8961554 ;
0086 13.52718 0.4809653 ;
0087 13.95531 0.8901121 ;
0088 14.56318 0.6767258 ;
0089 14.59449 1.190709 ;
0090 13.11096 0.2861982 ];
0091
0092 v_v=reshape([0.196497 -0.002605404 -0.002605404 0.05495016 ;
0093 0.6054919 0.007776652 0.007776652 0.01899244 ;
0094 0.5944617 0.0485788 0.0485788 0.03511229 ;
0095 0.3871268 0.0292966 0.0292966 0.02046839 ;
0096 0.3377683 0.02839657 0.02839657 0.04756354 ;
0097 1.00439 0.03595795 0.03595795 0.006737475 ]',[2 2 6]);
0098 end
0099
0100
0101
0102
0103
0104
0105 dpwtdef=[1.0000, 0.8250, 0.01868, 0.006773, 98.9, -0.4238];
0106
0107
0108 end
0109
0110
0111 p.fstep=5;
0112 p.fmax=4000;
0113 p.fres = 20;
0114 p.fbanklo = 10;
0115 p.mpsmooth = 21;
0116
0117 p.shortut = 7;
0118 p.pefact = 1.8;
0119 p.numopt = 3;
0120 p.flim = [60 400];
0121 p.w = dpwtdef;
0122
0123
0124 p.tmf = 2;
0125 p.tinc = 0.01;
0126 p.sopt = 'ilcwpf';
0127
0128
0129
0130 if nargin>=5 && isstruct(pp)
0131 fnq=fieldnames(pp);
0132 for i=1:length(fnq)
0133 if isfield(p,fnq{i})
0134 p.(fnq{i})=pp.(fnq{i});
0135 end
0136 end
0137 end
0138
0139
0140 if nargin>=3 && numel(tinc)>0
0141 p.tinc = tinc;
0142 end
0143 if nargin<4
0144 m='';
0145 end
0146
0147
0148 fmin = 0; fstep = p.fstep; fmax = p.fmax;
0149 fres = p.fres;
0150 [tx,f,MIX]=v_spgrambw(s,fs,fres,[fmin fstep fmax],[],p.tinc);
0151 nframes=length(tx);
0152 txinc=tx(2)-tx(1);
0153
0154
0155 [trans,cf]=v_filtbankm(length(f),2*length(f)-1,2*f(end),p.fbanklo,f(end),'usl');
0156 O = MIX*trans';
0157
0158
0159
0160
0161
0162 ltass = v_stdspectrum(6,'p',cf);
0163 auxf = [cf(1),(cf(1:end-1)+cf(2:end))./2,cf(end)];
0164 ltass = ltass.*diff(auxf);
0165
0166
0167 O = O.*repmat(diff(auxf),nframes,1);
0168 O1 = O;
0169
0170 if tx(end)<p.shortut
0171 eltass = mean(O,1);
0172 eltass = smooth(eltass,p.mpsmooth);
0173 eltass= eltass(:).';
0174
0175
0176 alpha = (ltass)./(eltass);
0177 alpha = alpha(:).';
0178 alpha = repmat(alpha,nframes,1);
0179 O = O.*alpha;
0180
0181
0182 else
0183
0184 tsmo = 3;
0185 stt = round(tsmo/txinc);
0186 eltass = timesm(O,stt);
0187 eltass = smooth(eltass,p.mpsmooth);
0188
0189
0190 alpha = repmat(ltass,nframes,1)./(eltass);
0191 O = O.*alpha;
0192
0193 end
0194
0195
0196
0197 ini = find(cf>3*cf(1));
0198 sca = cf/cf(ini(1));
0199
0200
0201 sca = sca(sca<10.5 & sca>0.5);
0202
0203 sca1 = sca;
0204 filh = 1./(p.pefact-cos(2*pi*sca1));
0205 filh = filh - sum(filh(1:end).*diff([sca1(1),(sca1(1:end-1)+sca1(2:end))./2,sca1(end)]))/sum(diff([sca1(1),(sca1(1:end-1)+sca1(2:end))./2,sca1(end)]));
0206
0207 posit = find(sca>=1);
0208 negat = find(sca<1);
0209 numz = length(posit)-1-length(negat);
0210 filh = filh./max(filh);
0211 filh = [zeros(1,numz) filh];
0212
0213
0214
0215 B = imfilter(O,filh);
0216
0217
0218 numopt = p.numopt;
0219 flim = p.flim;
0220 pfreq = find(cf>flim(1) & cf<flim(2));
0221 ff = zeros(nframes,numopt);
0222 amp = zeros(nframes,numopt);
0223 for i=1:nframes
0224 [pos,peak]=v_findpeaks(B(i,pfreq),[],5/(cf(pfreq(2))-cf(pfreq(1))));
0225 if numel(pos)
0226 [peak,ind]=sort(peak,'descend');
0227 pos = pos(ind);
0228 posff = cf(pfreq(pos));
0229 fin = min(numopt,length(posff));
0230 ff(i,1:fin)=posff(1:fin);
0231 amp(i,1:fin)=peak(1:fin);
0232 end
0233 end
0234
0235
0236
0237
0238
0239
0240
0241
0242 pow = mean(O,2);
0243
0244 vuvfea = [log(pow) 1e-3*sum(amp,2)./(pow+1.75*1e5)];
0245
0246
0247
0248 pru=v_gaussmixp(vuvfea,m_u,v_u,w_u);
0249 prv=v_gaussmixp(vuvfea,m_v,v_v,w_v);
0250
0251 pv=(1+exp(pru-prv)).^(-1);
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263 w = p.w;
0264
0265
0266 camp = -amp./repmat(max(amp,[],2),1,numopt);
0267 camp(amp==0)=w(5);
0268
0269
0270 tmf = p.tmf;
0271 inmf = round(tmf/txinc);
0272
0273
0274
0275
0276 cost = zeros(nframes,numopt);
0277 prev = zeros(nframes,numopt);
0278 medfx = zeros(nframes,1);
0279 dffact=2/txinc;
0280
0281
0282
0283 cost(1,:) = w(1)*camp(1,:);
0284 fpos = ff(1:min(inmf,end),1);
0285 mf=median(fpos(pv(1:min(inmf,end))>0.6));
0286 if isnan(mf)
0287 mf=median(fpos(pv(1:min(inmf,end))>0.5));
0288 if isnan(mf)
0289 mf=median(fpos(pv(1:min(inmf,end))>0.4));
0290 if isnan(mf)
0291 mf=median(fpos(pv(1:min(inmf,end))>0.3));
0292 if isnan(mf)
0293 mf=0;
0294 end
0295 end
0296 end
0297 end
0298 medfx(1)=mf;
0299
0300 for i=2:nframes
0301 if i>inmf
0302 fpos = ff(i-inmf:i,1);
0303 mf=median(fpos(pv(1:inmf)>0.6));
0304 if isnan(mf)
0305 mf=median(fpos(pv(1:inmf)>0.5));
0306 if isnan(mf)
0307 mf=median(fpos(pv(1:inmf)>0.4));
0308 if isnan(mf)
0309 mf=median(fpos(pv(1:inmf)>0.3));
0310 if isnan(mf)
0311 mf=0;
0312 end
0313 end
0314 end
0315 end
0316 end
0317 medfx(i)=mf;
0318
0319 df = dffact*(repmat(ff(i,:).',1,numopt) - repmat(ff(i-1,:),numopt,1))./(repmat(ff(i,:).',1,numopt) + repmat(ff(i-1,:),numopt,1));
0320 costdf=w(3)*min((df-w(6)).^2,w(4));
0321
0322
0323 if mf==0
0324 costf = zeros(1,numopt);
0325 else
0326 costf = abs(ff(i,:) - mf)./mf;
0327 end
0328 [cost(i,:),prev(i,:)]=min(costdf + repmat(cost(i-1,:),numopt,1),[],2);
0329 cost(i,:)=cost(i,:)+w(2)*costf + w(1)*camp(i,:);
0330
0331 end
0332
0333
0334
0335 fx=zeros(nframes,1);
0336 ax=zeros(nframes,1);
0337 best = zeros(nframes,1);
0338
0339 nose=find(cost(end,:)==min(cost(end,:)));
0340 best(end)=nose(1);
0341 fx(end)=ff(end,best(end));
0342 ax(end)=amp(end,best(end));
0343 for i=nframes:-1:2
0344 best(i-1)=prev(i,best(i));
0345 fx(i-1)=ff(i-1,best(i-1));
0346 ax(i-1)=amp(i-1,best(i-1));
0347 end
0348
0349 if nargout>=4
0350 fv.vuvfea=vuvfea;
0351 fv.best=best;
0352 fv.ff=ff;
0353 fv.amp=amp;
0354 fv.medfx=medfx;
0355 fv.w=w;
0356 fv.dffact=dffact;
0357 fv.hist = [log(mean(O,2)) sum(amp,2)./((mean(O,2)))];
0358 end
0359
0360 if ~nargout || any(m=='g') || any(m=='G')
0361 nax=0;
0362 msk=pv>0.5;
0363 fxg=fx;
0364 fxg(~msk)=NaN;
0365 fxb=fx;
0366 fxb(msk)=NaN;
0367 if any(m=='G') || ~nargout && ~any(m=='g')
0368 clf;
0369 v_spgrambw(s,fs,p.sopt);
0370 hold on
0371 plot(tx,log10(fxg),'-b',tx,log10(fxb),'-r');
0372 yy=get(gca,'ylim');
0373 plot(tx,yy(1)+yy*[-1;1]*(0.02+0.05*pv),'-k');
0374 hold off
0375 nax=nax+1;
0376 axh(nax)=gca;
0377 if any(m=='g')
0378 figure;
0379 end
0380 end
0381 if any(m=='g')
0382 ns=length(s);
0383 [tsr,ix]=sort([(1:ns)/fs 0.5*(tx(1:end-1)+tx(2:end))']);
0384 jx(ix)=1:length(ix);
0385 sp2fr=jx(1:ns)-(0:ns-1);
0386 spmsk=msk(sp2fr);
0387 sg=s;
0388 sg(~spmsk)=NaN;
0389 sb=s;
0390 sb(spmsk)=NaN;
0391 clf;
0392 subplot(5,1,1);
0393 plot(tx,pv,'-b',(1:ns)/fs,0.5*mod(cumsum(fx(sp2fr)/fs),1)-0.6,'-b');
0394 nax=nax+1;
0395 axh(nax)=gca;
0396 ylabel('\phi(t), P(V)');
0397 set(gca,'ylim',[-0.65 1.05]);
0398 subplot(5,1,2:3);
0399 plot((1:ns)/fs,sg,'-b',(1:ns)/fs,sb,'-r');
0400 nax=nax+1;
0401 axh(nax)=gca;
0402 subplot(5,1,4:5);
0403 plot(tx,fxg,'-b',tx,fxb,'-r');
0404 ylabel('Pitch (Hz)');
0405
0406
0407 set(gca,'ylim',[min(fxg)-30 max(fxg)+30]);
0408 nax=nax+1;
0409 axh(nax)=gca;
0410 end
0411 if nax>1
0412 linkaxes(axh,'x');
0413 end
0414 end
0415
0416 function y=smooth(x,n)
0417 nx=size(x,2);
0418 nf=size(x,1);
0419 c=cumsum(x,2);
0420 y=[c(:,1:2:n)./repmat(1:2:n,nf,1) (c(:,n+1:end)-c(:,1:end-n))/n (repmat(c(:,end),1,floor(n/2))-c(:,end-n+2:2:end-1))./repmat(n-2:-2:1,nf,1)];
0421
0422 function y=timesm(x,n)
0423 if ~mod(n,2)
0424 n = n+1;
0425 end
0426 nx=size(x,2);
0427 nf=size(x,1);
0428 c=cumsum(x,1);
0429 mid = round(n/2);
0430 y=[c(mid:n,:)./repmat((mid:n).',1,nx); ...
0431 (c(n+1:end,:)-c(1:end-n,:))/n; ...
0432 (repmat(c(end,:),mid-1,1) - c(end-n+1:end-mid,:))./repmat((n-1:-1:mid).',1,nx)];