0001 function [m,v]=v_psycdigit(proc,r,mode,p,q,xp,noise,fn,dfile,ofile)
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
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121 persistent tok mtok tigflgp digitsp
0122 if nargin<10
0123 ofile='';
0124 if nargin<9
0125 dfile='F:\home\dmb\data\old_speech\tidigits';
0126
0127 if nargin<8
0128 fn=16000;
0129 if nargin<7
0130 noise=[];
0131 if nargin<6
0132 xp=[];
0133 if nargin<5
0134 q=[];
0135 if nargin<4
0136 p=[];
0137 if nargin<3
0138 mode='';
0139 if nargin<2
0140 r=[];
0141 end
0142 end
0143 end
0144 end
0145 end
0146 end
0147 end
0148 end
0149 end
0150
0151
0152
0153 i=1;
0154 pv=[0 3 3 1 25 0 0 1 1 1 1 100];
0155 px=[0 3 3 1 25 20 18 1 1 1 1 100];
0156 pvs='jlLmnxXvVfFb';
0157 pf=zeros(1,52);
0158 mval=0;
0159 lmode=length(mode);
0160 while i<=lmode
0161 if i<lmode
0162 [v,nv,e,ni]=sscanf(mode(i+1:end),'%d',1);
0163 else
0164 nv=0;
0165 ni=1;
0166 end
0167 j=find(mode(i)==pvs,1);
0168 k=1+2*(double(lower(mode(i)))-'a')+(mode(i)<'a');
0169 if k>=1 && k<=52
0170 pf(k)=1-pf(k);
0171 if ~isempty(j)
0172 if nv==0
0173 pv(j)=px(j);
0174 else
0175 pv(j)=v;
0176 mval=mval || j==4;
0177 end
0178 end
0179 end
0180 i=i+ni;
0181 end
0182 if isempty(proc)
0183 pv(4)=0;
0184 end
0185
0186
0187 varthr=20;
0188 nxevo=30;
0189
0190 if pf(23)>pf(24) || pv(2)>pv(3)
0191 pv(3)=pv(2);
0192 end
0193 if pf(24)>pf(23)
0194 pv(2)=1;
0195 end
0196 zmodel=pv(4)+(pf(26)>0);
0197 ntrial=zmodel*pv(5);
0198 if pf(12)>pf(11)
0199 pv(10)=pv(11);
0200 end
0201 if pf(11)+pf(12)>0 && (pv(10)<1 || pv(10)>3)
0202 error('Invalid average type for option f/F');
0203 end
0204
0205
0206 unimode=pf(41);
0207 if isempty(p)
0208 p=0.75;
0209 end
0210 if size(p,1)<3-unimode
0211 p(2-unimode,1)=0.04;
0212 p(3-unimode,1)=0;
0213 end
0214 if size(p,2)<(zmodel)
0215 p=p(:,1+mod(0:zmodel-1,size(p,2)));
0216 end
0217
0218
0219
0220 tigflg=[pv(2:3) pf(51:52)];
0221 if isempty(tok) || any(tigflg~=tigflgp) || ~strcmp(dfile,digitsp)
0222 disp('Scanning TIDIGITS files ...');
0223 digitsp=dfile;
0224 tigflgp=tigflg;
0225 if any(dfile(end)=='/\')
0226 dfile(end)=[];
0227 end
0228 dirlist{1}='';
0229 ntok=0;
0230 mtok=zeros(1,pv(3));
0231 tok=cell(1,2);
0232 while ~isempty(dirlist)
0233 dd=dir([dfile dirlist{1}]);
0234 for i=1:length(dd)
0235 name=dd(i).name;
0236 if name(1)~='.'
0237 if dd(i).isdir
0238 dirlist{end+1}=[dirlist{1} '\' name];
0239 elseif length(name)>4 && strcmpi(name(end-3:end),'.wav')
0240 digs=name(1:end-4);
0241 digz=upper(digs)=='Z';
0242 digo=upper(digs)=='O';
0243 digs(digo | digz)='0';
0244 digs=digs(digs>='0' & digs<='9');
0245 ndigs=length(digs);
0246 if ndigs>=tigflg(1) && ndigs<=tigflg(2) && ~(tigflg(3) && any(digo)) && ~(tigflg(4) && any(digz))
0247 ntok=ntok+1;
0248 mtok(ndigs)=mtok(ndigs)+1;
0249 tok{ntok,1}=[dirlist{1} '\' name];
0250 tok{ntok,2}=digs;
0251 end
0252 end
0253 end
0254 end
0255 dirlist(1)=[];
0256 end
0257 end
0258 ntok=size(tok,1);
0259 if pf(46)
0260 [s,fs]=v_readwav([dfile tok{1,1}]);
0261 else
0262 [s,fs]=v_readsph([dfile tok{1,1}]);
0263 end
0264
0265
0266
0267 if any(mode=='d')
0268 p(3-unimode,:)=0.1;
0269 else
0270 p(3-unimode,:)=0.1.^(1:pv(3))*mtok'/ntok;
0271 end
0272
0273
0274
0275 if unimode
0276 [xx,ii,m,v]=v_psycestu(-zmodel,p,q,xp);
0277 [pact,qact]=v_psycestu(0);
0278 else
0279 [xx,ii,m,v]=v_psycest(-zmodel,p,q,xp);
0280 [pact,qact]=v_psycest(0);
0281 end
0282 if pv(4) && pf(5)
0283 x=[];
0284 ii=1;
0285 if pf(1)
0286 if isempty(r)
0287 if mval
0288 procdesc=proc(x,fs,xx,ii);
0289 else
0290 procdesc=proc(x,fs,xx);
0291 end
0292 else
0293 if mval
0294 procdesc=proc(x,fs,xx,ii,r);
0295 else
0296 procdesc=proc(x,fs,xx,r);
0297 end
0298 end
0299 else
0300 if isempty(r)
0301 if mval
0302 procdesc=proc(x,fs,ii);
0303 else
0304 procdesc=proc(x,fs);
0305 end
0306 else
0307 if mval
0308 procdesc=proc(x,fs,ii,r);
0309 else
0310 procdesc=proc(x,fs,r);
0311 end
0312 end
0313 end
0314 else
0315 procdesc='';
0316 end
0317
0318
0319
0320 if pf(29) || pf(30)
0321 nw=fix(datevec(now));
0322 of=sprintf('psy%4d%02d%02d%02d%02d.txt',nw(1:5));
0323 ofid=fopen(of,'wt');
0324 if ~ofid
0325 error('Cannot write to %s',of);
0326 end
0327 fprintf(ofid,'%% %s evaluation on %s\n',mfilename,datestr(now));
0328 fmfnm=[mfilename('fullpath') '.m'];
0329 dd=dir(fmfnm);
0330 fprintf(ofid,'%% %s = %d bytes = %s\n',fmfnm,dd.bytes,dd.date);
0331 fprintf(ofid,'V %d\n',2);
0332 fmfnm=which(func2str(proc));
0333 dd=dir(fmfnm);
0334 fprintf(ofid,'P %s = %d bytes = %s\n',fmfnm,dd.bytes,dd.date);
0335 if numel(procdesc)
0336 fprintf(ofid,'C %s\n',procdesc);
0337 end
0338 fprintf(ofid,'O %s\n',mode);
0339 end
0340
0341
0342
0343 disp(['Testing ' procdesc]);
0344
0345 mnt=zeros(2+2*unimode,zmodel,3,ntrial+1);
0346 vnt=zeros(3+unimode,zmodel,ntrial+1);
0347 mnt(:,:,:,1)=m;
0348 vnt(:,:,1)=v;
0349 i=0;
0350 imax=0;
0351 quitit=0;
0352 while ~quitit
0353 i=i+1;
0354 isp=min(1+floor(rand(1)*ntok),ntok);
0355 if pf(46)
0356 [s,fs]=v_readwav([dfile tok{isp,1}]);
0357 else
0358 [s,fs]=v_readsph([dfile tok{isp,1}]);
0359 end
0360 s=[zeros(pv(6)*round(fs/10),1); v_activlev(s(:),fs,'n')];
0361 if pf(1) && ii<=pv(4)
0362 if isempty(r)
0363 if mval
0364 y=proc(s,fs,xx,ii);
0365 else
0366 y=proc(s,fs,xx);
0367 end
0368 else
0369 if mval
0370 y=proc(s,fs,xx,ii,r);
0371 else
0372 y=proc(s,fs,xx,r);
0373 end
0374 end
0375 else
0376 nn=randn(length(s),1);
0377 x=nn+s*10^(xx/20);
0378 if ii<=pv(4)
0379 if isempty(r)
0380 if mval
0381 y=proc(x,fs,ii);
0382 else
0383 y=proc(x,fs);
0384 end
0385 else
0386 if mval
0387 y=proc(x,fs,ii,r);
0388 else
0389 y=proc(x,fs,r);
0390 end
0391 end
0392 else
0393 y=x;
0394 end
0395 end
0396 y=y(1+pv(7)*round(fs/10):end);
0397 if pf(13)
0398 prg=sprintf(' %d',length(tok{isp,2}));
0399 else
0400 prg='';
0401 end
0402 if pf(14)
0403 prG=sprintf('SNR=%d dB, ',round(xx));
0404 else
0405 prG='';
0406 end
0407 if pf(35)
0408 prr=sprintf(', ''r'' to repeat');
0409 else
0410 prr='';
0411 end
0412 prompt=[prG 'enter' prg ' digits (''q'' to quit' prr '): '];
0413
0414 ansr=-1;
0415 say=1;
0416 while ansr<0
0417 if say
0418 tdel=0;
0419 tic;
0420 soundsc(y,fs);
0421 say=0;
0422 end
0423 rv=input(prompt,'s');
0424 tdel=toc;
0425 if ~isempty(rv)
0426 if lower(rv(1))=='q'
0427 quitit=1;
0428 ansr=2;
0429 elseif lower(rv(1))=='r' && pf(35)
0430 say=1;
0431 elseif lower(rv(1))=='s' && pf(37)
0432 ofn=input('File name: ','s');
0433 if numel(ofn)
0434 v_writewav(y,fs,ofn);
0435 end
0436 elseif all(rv>='0') && all(rv<='9') && ( ~pf(13) || length(rv)==length(tok{isp,2}))
0437 ansr=strcmp(rv,tok{isp,2});
0438 end
0439 end
0440 end
0441 quitit=quitit || i==ntrial;
0442 jj=ii;
0443 xxold=xx;
0444 if ansr<2
0445 if unimode
0446 [xx,ii,m,v]=v_psycestu(ii,xx,ansr);
0447 else
0448 [xx,ii,m,v]=v_psycest(ii,xx,ansr);
0449 end
0450 mnt(:,:,:,i+1)=m;
0451 vnt(:,:,i+1)=v;
0452 imax=i;
0453 end
0454 if pf(30) || quitit && pf(29)
0455 if ansr>1
0456 rv=num2str(double(rv(1)));
0457 end
0458
0459 fprintf(ofid,'M %d %d %.3g %s %s %d %.1f',i,ii,xxold,tok{isp,2},rv,ansr,tdel);
0460 fprintf(ofid,' %.3g',m(:,jj,:),v(:,jj));
0461 fprintf(ofid,'\n');
0462 end
0463 if pf(32)
0464 figure(jj);
0465 if unimode
0466 v_psycestu(jj);
0467 else
0468 v_psycest(jj);
0469 end
0470 elseif quitit && pf(31)
0471 for jj=1:zmodel
0472 figure(jj);
0473 if unimode
0474 v_psycestu(jj);
0475 else
0476 v_psycest(jj);
0477 end
0478 end
0479 end
0480 if pf(12) || quitit && pf(11)
0481 figure(pv(12)+1);
0482 if unimode
0483 qqu.pk=m(1,jj,pv(10));
0484 qqu.w=m(2,1,pv(10));
0485 qqu.ll=m(3,1,pv(10));
0486 qqu.lh=m(4,1,pv(10));
0487 qqu.gu=pact(2);
0488 qqu.la=pact(1);
0489 sw=qqu.w*qqu.ll*qqu.lh/(qqu.ll+qqu.lh);
0490 xax=linspace(qqu.pk-(4+sw)/qqu.ll,qqu.pk+(4+sw)/qqu.lh,200);
0491 bs=(qqu.lh-qqu.ll)/2;
0492 bu=(qqu.lh+qqu.ll)/2;
0493 plot(xax,qqu.gu+(1-qqu.gu-qqu.la)*(1+2*exp(-sw)*cosh(bs*(xax-qqu.pk)+bu*abs(xax-qqu.pk))+exp(-2*sw)).^(-1));
0494 else
0495 sd=(pact(1,:)-pact(3,:)).*(1-pact(2,:)-pact(1,:))./(m(2,:,pv(10)).*(1-pact(3,:)-pact(2,:)));
0496 md=m(1,:,pv(10))-sd.*log((pact(1,:)-pact(3,:))./(1-pact(2,:)-pact(1,:)));
0497 xax=linspace(min(md-3*sd),max(md+3*sd),100);
0498 for jj=1:zmodel
0499 plot(xax,v_psychofunc('',[pact(1,jj); m(:,jj,pv(10)); pact(2:3,jj); qact(10)],xax));
0500 hold on
0501 end
0502 hold off
0503 end
0504 axis([xax(1) xax(end) 0 1]);
0505 xlabel('SNR (dB)');
0506 ylabel('Recognition Probability');
0507 title(sprintf('Intelligibility: %s',procdesc));
0508 end
0509 if pf(10) || quitit && pf(9)
0510 figure(pv(12)+3);
0511 psyevo=zeros(ntrial+1,nxevo);
0512 if unimode
0513 qqu.pk=m(1,jj,pv(10));
0514 qqu.w=m(2,1,pv(10));
0515 qqu.ll=m(3,1,pv(10));
0516 qqu.lh=m(4,1,pv(10));
0517 qqu.gu=pact(2);
0518 qqu.la=pact(1);
0519 sw=qqu.w*qqu.ll*qqu.lh/(qqu.ll+qqu.lh);
0520 xax=linspace(qqu.pk-(4+sw)/qqu.ll,qqu.pk+(4+sw)/qqu.lh,nxevo);
0521 for iet=1:imax+1
0522 qqu.pk=mnt(1,1,pv(10),iet);
0523 qqu.w=mnt(2,1,pv(10),iet);
0524 qqu.ll=mnt(3,1,pv(10),iet);
0525 qqu.lh=mnt(4,1,pv(10),iet);
0526 sw=qqu.w*qqu.ll*qqu.lh/(qqu.ll+qqu.lh);
0527 bs=(qqu.lh-qqu.ll)/2;
0528 bu=(qqu.lh+qqu.ll)/2;
0529 psyevo(iet,:)=qqu.gu+(1-qqu.gu-qqu.la)*(1+2*exp(-sw)*cosh(bs*(xax-qqu.pk)+bu*abs(xax-qqu.pk))+exp(-2*sw)).^(-1);
0530 end
0531 imagesc(xax,0:ntrial,psyevo);
0532 axis 'ij'
0533 else
0534 end
0535 xlabel('SNR (dB)');
0536 ylabel('After trial');
0537 title(sprintf('Intelligibility: %s',procdesc));
0538 end
0539 if pf(44) || quitit && pf(43)
0540 figure(pv(12)+2);
0541 if unimode
0542 sw=mnt(2,1,pv(10),1:imax+1).*mnt(3,1,pv(10),1:imax+1).*mnt(4,1,pv(10),1:imax+1)./(mnt(3,1,pv(10),1:imax+1).*mnt(4,1,pv(10),1:imax+1));
0543 subplot(221)
0544 plot(0:imax,squeeze(mnt(1,1,pv(10),1:imax+1)));
0545 set(gca,'xlim',[0 ntrial]);
0546 xlabel('After trial');
0547 ylabel('Peak position (dB SNR)');
0548 subplot(222)
0549 plot(0:imax,squeeze(mnt(1,1,pv(10),1:imax+1)-sw.*mnt(3,1,pv(10),1:imax+1)),'-b',0:imax,squeeze(mnt(1,1,pv(10),1:imax+1)+sw.*mnt(4,1,pv(10),1:imax+1)),'-b');
0550 set(gca,'xlim',[0 ntrial]);
0551 xlabel('After trial');
0552 ylabel('Inflections (dB)');
0553 subplot(223)
0554 plot(0:imax,squeeze(mnt(3,1,pv(10),1:imax+1)));
0555 set(gca,'xlim',[0 ntrial]);
0556 xlabel('After trial');
0557 ylabel('Upwards slope (prob/dB)');
0558 subplot(224)
0559 plot(0:imax,squeeze(mnt(4,1,pv(10),1:imax+1)));
0560 set(gca,'xlim',[0 ntrial]);
0561 xlabel('After trial');
0562 ylabel('Downwards slope (prob/dB)');
0563 else
0564 subplot(211);
0565 for jj=1:zmodel
0566 plot(0:imax,squeeze(mnt(1,jj,pv(10),1:imax+1)));
0567 hold on
0568 end
0569 hold off
0570 set(gca,'xlim',[0 ntrial]);
0571 xlabel('After trial');
0572 ylabel('SRT (dB SNR)');
0573 subplot(212);
0574 for jj=1:zmodel
0575 plot(0:i,squeeze(mnt(2,jj,pv(10),1:imax+1)));
0576 hold on
0577 end
0578 hold off
0579 set(gca,'xlim',[0 ntrial]);
0580 xlabel('After trial');
0581 ylabel('Slope (Prob/dB)');
0582 end
0583 end
0584 end
0585 if pf(29) || pf(30)
0586 fclose(ofid);
0587 end