Home > voicebox > psycdigit.m

psycdigit

PURPOSE ^

PSYCDIGIT measures psychometric function using TIDIGITS stimuli

SYNOPSIS ^

function [m,v]=psycdigit(proc,r,mode,p,q,xp,noise,fn,dfile,ofile)

DESCRIPTION ^

PSYCDIGIT measures psychometric function using TIDIGITS stimuli

 Usage:
         (1)[m,v]=psycdigit([],[],'GMPWrn10',[],[],[],[],[],dfile);
                       % measure SRT using addditive white noise, repetitions allowed, data in .WAV format
         (2)[m,v]=psycdigit(@specsub,[],'GMPn10',[],[],[],[],[],dfile);
                       % compare spectral subtraction with unprocessed noisy speech

 Inputs:
         proc    processing function handle (e.g. @specsub) called as follows:
                    y=proc(x,fs,i,r)  % process noisy waveform x through model i with parameters r
                    y=proc(x,fs,snr,i,r) % process clean speech degraded to snr through model i with parameters r
                    y=proc()   % return comment text string describing algorithm (if 'c' option specified)
         r       parameters for proc (omitted from call if r=[])
         mode    string containing options (see below for list)
         p,q,xp  parameters passed on to psycest or psycestu
         noise   noise waveform or wav file containing noise
         fn      noise waveform sample frequency [16000]
         dfile   path to TIdigits folder
         ofile   output text file (see below for output file format)

 Outputs:
          m(2,n,3) estimated srt and slope of all models
                   m(i,n,k): i={1=srt, 2=slope}, n=model number, k={1=mean, 2=mode (MAP), 3=marginal mode}
          v(3,n)   estimated covariance matrix entries:
                   [var(srt) cov(srt,slope) var(slope)]' of n'th model

 List of options for "mode" input:
         a       proc adds its own noise
         b*      base figure number for plotting results [100]
         c       y=proc() returns comment string
         e/E*    plot evolving psychometric functions *=1,2,3 for mean, mode, marginal mode (F=after each trial)
         f/F*    plot psychometric functions *=1,2,3 for mean, mode, marginal mode (F=after each trial)
         g       prompt with number of digits
         G       prompt with SNR
         l*      min token length (in digits)
         L*      max token length (if different from min)
         m*      use * external models [default 1]
         M       include an extra model with no processing
         n*      *=average number of trials per model
         p/P     plot pdf (P=after each trial)
         r       allow repetitions
         s       respond s to save the noisy stimulus as a .wav file
         t/T*    taper noise */10 seconds at start and end
         v/V*    plot srt/slope convergence (V= after each trial)
                 *=1,2,3 for mean, mode, marginal mode
         x*      add */10 seconds of noise to the front of the speech before processing
         X*      truncate */10 seconds from the start of the processed speech
         W       data is in microsoft WAV format
         z       omit tokens containing "oh"
         Z       omit tokens containing "zero"

 Output file format
   Each line starts 'x ' where x is one of the following
       %  Comment
       V  File version type
       O  mode options
       P  details about proc
       C  Comment returned by proc
       M  measurement

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [m,v]=psycdigit(proc,r,mode,p,q,xp,noise,fn,dfile,ofile)
0002 %PSYCDIGIT measures psychometric function using TIDIGITS stimuli
0003 %
0004 % Usage:
0005 %         (1)[m,v]=psycdigit([],[],'GMPWrn10',[],[],[],[],[],dfile);
0006 %                       % measure SRT using addditive white noise, repetitions allowed, data in .WAV format
0007 %         (2)[m,v]=psycdigit(@specsub,[],'GMPn10',[],[],[],[],[],dfile);
0008 %                       % compare spectral subtraction with unprocessed noisy speech
0009 %
0010 % Inputs:
0011 %         proc    processing function handle (e.g. @specsub) called as follows:
0012 %                    y=proc(x,fs,i,r)  % process noisy waveform x through model i with parameters r
0013 %                    y=proc(x,fs,snr,i,r) % process clean speech degraded to snr through model i with parameters r
0014 %                    y=proc()   % return comment text string describing algorithm (if 'c' option specified)
0015 %         r       parameters for proc (omitted from call if r=[])
0016 %         mode    string containing options (see below for list)
0017 %         p,q,xp  parameters passed on to psycest or psycestu
0018 %         noise   noise waveform or wav file containing noise
0019 %         fn      noise waveform sample frequency [16000]
0020 %         dfile   path to TIdigits folder
0021 %         ofile   output text file (see below for output file format)
0022 %
0023 % Outputs:
0024 %          m(2,n,3) estimated srt and slope of all models
0025 %                   m(i,n,k): i={1=srt, 2=slope}, n=model number, k={1=mean, 2=mode (MAP), 3=marginal mode}
0026 %          v(3,n)   estimated covariance matrix entries:
0027 %                   [var(srt) cov(srt,slope) var(slope)]' of n'th model
0028 %
0029 % List of options for "mode" input:
0030 %         a       proc adds its own noise
0031 %         b*      base figure number for plotting results [100]
0032 %         c       y=proc() returns comment string
0033 %         e/E*    plot evolving psychometric functions *=1,2,3 for mean, mode, marginal mode (F=after each trial)
0034 %         f/F*    plot psychometric functions *=1,2,3 for mean, mode, marginal mode (F=after each trial)
0035 %         g       prompt with number of digits
0036 %         G       prompt with SNR
0037 %         l*      min token length (in digits)
0038 %         L*      max token length (if different from min)
0039 %         m*      use * external models [default 1]
0040 %         M       include an extra model with no processing
0041 %         n*      *=average number of trials per model
0042 %         p/P     plot pdf (P=after each trial)
0043 %         r       allow repetitions
0044 %         s       respond s to save the noisy stimulus as a .wav file
0045 %         t/T*    taper noise */10 seconds at start and end
0046 %         v/V*    plot srt/slope convergence (V= after each trial)
0047 %                 *=1,2,3 for mean, mode, marginal mode
0048 %         x*      add */10 seconds of noise to the front of the speech before processing
0049 %         X*      truncate */10 seconds from the start of the processed speech
0050 %         W       data is in microsoft WAV format
0051 %         z       omit tokens containing "oh"
0052 %         Z       omit tokens containing "zero"
0053 %
0054 % Output file format
0055 %   Each line starts 'x ' where x is one of the following
0056 %       %  Comment
0057 %       V  File version type
0058 %       O  mode options
0059 %       P  details about proc
0060 %       C  Comment returned by proc
0061 %       M  measurement
0062 
0063 % Future mode options:
0064 %        [ d     score as single digits ]
0065 %        [ i/I   plot SRT improvement (I=after each trial) ]
0066 %        [ j*    scaling: 0=autoscale each token, 1=constant speech,2=const noise, 3=const noisy ]
0067 %        [N*     ignore the first * trials ]
0068 %        [o/O    write to output file, O write result of each probe]
0069 %        [ S     save all stimuli as wav files ]
0070 %        [ u     do not normalize the noise level ]
0071 %        [ y*    type of noise ]
0072 %
0073 % Bugs/Suggestions
0074 % (1) Add sounds to indicate error and end of test
0075 % (2) Routine could return a label to replace "SNR" in necessary
0076 % (3) Add an input option to "discard" this sample
0077 % (4) output file should include mode argument + date + IP address + computer name + r arguments + p/q values
0078 % (5) Quote filenames in output file or else never have anything else on the line
0079 % (6) silence detection doesn't work
0080 
0081 %      Copyright (C) Mike Brookes 2010
0082 %      Version: $Id: psycdigit.m 8166 2016-07-08 14:46:15Z dmb $
0083 %
0084 %   VOICEBOX is a MATLAB toolbox for speech processing.
0085 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
0086 %
0087 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0088 %   This program is free software; you can redistribute it and/or modify
0089 %   it under the terms of the GNU General Public License as published by
0090 %   the Free Software Foundation; either version 2 of the License, or
0091 %   (at your option) any later version.
0092 %
0093 %   This program is distributed in the hope that it will be useful,
0094 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0095 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0096 %   GNU General Public License for more details.
0097 %
0098 %   You can obtain a copy of the GNU General Public License from
0099 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0100 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
0101 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0102 
0103 %     25+4     m*      use * external models [default 1]
0104 %     26    M       include an extra model with no processing
0105 %      1   a       proc adds its own noise
0106 %      27+5   n*      *=average number of trials per model
0107 %     23+2    l*      min token length (in digits)
0108 %      24+3   L*      max token length (if different from min)
0109 %      7   d       score as single digits
0110 %      37   s       speech-shaped noise instead of white
0111 %       19+1  j*      scaling: 0=autoscale each token, 1=constant speech, 2=const noise, 3=const noisy
0112 %      35   r       allow repetitions
0113 %     41    u       unimodal psychometric function
0114 %     32    P       plot pdf
0115 %     31   p       plot srt convergence
0116 %     13    g       prompt with number of digits
0117 %      14   G       prompt with SNR
0118 %      47+6   x*      add */10 seconds of noise to the front of the speech
0119 %      48+7   X*      truncate */10 from the start of the processed speech
0120 % input parameters
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         %         dfile='Y:\Speech\TIDigits\Disc 1\TIDIGITS';
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 % parse the mode string
0152 
0153 i=1;
0154 pv=[0 3 3 1 25  0  0 1 1 1 1 100]; % pv default values if option is unspecified
0155 px=[0 3 3 1 25 20 18 1 1 1 1 100]; % pv values if option is given without a value
0156 pvs='jlLmnxXvVfFb';
0157 pf=zeros(1,52);
0158 mval=0;
0159 lmode=length(mode);
0160 while i<=lmode
0161     if i<lmode  % read a following integer if it exists
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;  % m has a value specified
0177             end
0178         end
0179     end
0180     i=i+ni;
0181 end
0182 if isempty(proc)
0183     pv(4)=0; % not allowed any processed models if not process is specified
0184 end
0185 % derived input parameters
0186 
0187 varthr=20;   % variance threshold for "silence"
0188 nxevo=30; % size of evloving pdf
0189 
0190 if pf(23)>pf(24) || pv(2)>pv(3)
0191     pv(3)=pv(2);                % Make max token length reasonable
0192 end
0193 if pf(24)>pf(23)
0194     pv(2)=1;                    % min =1 if only max is specified
0195 end
0196 zmodel=pv(4)+(pf(26)>0); % total number of models (including null model)
0197 ntrial=zmodel*pv(5);  % total number of trials
0198 if pf(12)>pf(11)    % 'F' specified but not 'f'
0199     pv(10)=pv(11);  % copy across the average type selector
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 % sort out p argument to psycest or psycestu
0206 unimode=pf(41);
0207 if isempty(p)
0208     p=0.75;     % default recognition rate at threshold
0209 end
0210 if size(p,1)<3-unimode
0211     p(2-unimode,1)=0.04;     % miss probability
0212     p(3-unimode,1)=0;     % dummy entry for guess probability
0213 end
0214 if size(p,2)<(zmodel)
0215     p=p(:,1+mod(0:zmodel-1,size(p,2)));
0216 end
0217 
0218 % first sort out the digit samples
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)=[];  % remove a trailing separator
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)~='.'   % ignore directories starting with '.'
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)=[];  % remove this directory from the list
0256     end
0257 end
0258 ntok=size(tok,1);
0259 if pf(46)
0260     [s,fs]=readwav([dfile tok{1,1}]); % get the first speech token to set fs
0261 else
0262     [s,fs]=readsph([dfile tok{1,1}]); % get the first speech token to set fs
0263 end
0264 
0265 % calculate guess probability assuming you get the correct number of digits
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 % now initialize the models
0274 
0275 if unimode
0276     [xx,ii,m,v]=psycestu(-zmodel,p,q,xp); % initialize all models
0277     [pact,qact]=psycestu(0);  % save the actual parameters
0278 else
0279     [xx,ii,m,v]=psycest(-zmodel,p,q,xp); % initialize all models
0280     [pact,qact]=psycest(0);  % save the actual parameters
0281 end
0282 if pv(4) && pf(5)
0283     x=[];  % set x=[] to force a description output
0284     ii=1;   % for now, only do model 1
0285     if pf(1)   % if process adds its own noise
0286         if isempty(r)  % proc does not require any auxilliary parameters
0287             if mval
0288                 procdesc=proc(x,fs,xx,ii);  % process the noisy speech
0289             else
0290                 procdesc=proc(x,fs,xx);  % process the noisy speech
0291             end
0292         else
0293             if mval
0294                 procdesc=proc(x,fs,xx,ii,r);  % process the noisy speech
0295             else
0296                 procdesc=proc(x,fs,xx,r);  % process the noisy speech
0297             end
0298         end
0299     else   % if process does not add its own noise
0300         if isempty(r)  % proc does not require any auxilliary parameters
0301             if mval
0302                 procdesc=proc(x,fs,ii);  % process the noisy speech
0303             else
0304                 procdesc=proc(x,fs);  % process the noisy speech
0305             end
0306         else
0307             if mval
0308                 procdesc=proc(x,fs,ii,r);  % process the noisy speech
0309             else
0310                 procdesc=proc(x,fs,r);  % process the noisy speech
0311             end
0312         end
0313     end
0314 else
0315     procdesc='';
0316 end
0317 
0318 % now initialize the output file
0319 
0320 if pf(29) || pf(30)   % o/O output info
0321     nw=fix(datevec(now));
0322     of=sprintf('psy%4d%02d%02d%02d%02d.txt',nw(1:5)); % filename includes date and time to nearest minute
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); % print file format version number
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 % now start testing
0342 
0343 disp(['Testing ' procdesc]);
0344 % now do the main loop
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); % select a token
0355     if pf(46)
0356         [s,fs]=readwav([dfile tok{isp,1}]); % get the speech token
0357     else
0358         [s,fs]=readsph([dfile tok{isp,1}]); % get the speech token
0359     end
0360     s=[zeros(pv(6)*round(fs/10),1); activlev(s(:),fs,'n')]; % preappend zeros and normalize speech level
0361     if pf(1) && ii<=pv(4)  % if process adds its own noise
0362         if isempty(r)
0363             if mval
0364                 y=proc(s,fs,xx,ii);  % process the noisy speech
0365             else
0366                 y=proc(s,fs,xx);  % process the noisy speech
0367             end
0368         else
0369             if mval
0370                 y=proc(s,fs,xx,ii,r);  % process the noisy speech
0371             else
0372                 y=proc(s,fs,xx,r);  % process the noisy speech
0373             end
0374         end
0375     else
0376         nn=randn(length(s),1);
0377         x=nn+s*10^(xx/20);   % create the data
0378         if ii<=pv(4)
0379             if isempty(r)
0380                 if mval
0381                     y=proc(x,fs,ii);  % process the noisy speech
0382                 else
0383                     y=proc(x,fs);  % process the noisy speech
0384                 end
0385             else
0386                 if mval
0387                     y=proc(x,fs,ii,r);  % process the noisy speech
0388                 else
0389                     y=proc(x,fs,r);  % process the noisy speech
0390                 end
0391             end
0392         else
0393             y=x;            % unprocessed for last model
0394         end
0395     end
0396     y=y(1+pv(7)*round(fs/10):end);   % remove junk from the start
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     %     ansr=-(var(y)>varthr);  % meant to detect silences but doesn't work
0414     ansr=-1;
0415     say=1;
0416     while ansr<0
0417         if say
0418             tdel=0;
0419             tic;
0420             soundsc(y,fs);      % *** probably shouldn't be ...sc
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)   % save the token
0432                 ofn=input('File name: ','s');
0433                 if numel(ofn)
0434                     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;   % mark as quiting if we have done all the trials
0442     jj=ii;  % remember which model has just been updated
0443     xxold=xx; % and what the SNR was
0444     if ansr<2  % valid answer: update the pdfs
0445         if unimode
0446             [xx,ii,m,v]=psycestu(ii,xx,ansr);
0447         else
0448             [xx,ii,m,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)       % 'O/o': output
0455         if ansr>1
0456             rv=num2str(double(rv(1)));
0457         end
0458         % could add in token name and length
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)                           % 'P': plot PDF: figures 1:m
0464         figure(jj);
0465         if unimode
0466             psycestu(jj);
0467         else
0468             psycest(jj);
0469         end
0470     elseif quitit && pf(31)             % 'p': plot PDF: figures 1:m
0471         for jj=1:zmodel
0472             figure(jj);
0473             if unimode
0474                 psycestu(jj);
0475             else
0476                 psycest(jj);
0477             end
0478         end
0479     end
0480     if pf(12) || quitit && pf(11)       % 'F/f': plot Psychometric function on figure 101
0481         figure(pv(12)+1);
0482         if unimode
0483             qqu.pk=m(1,jj,pv(10));      % peak position
0484             qqu.w=m(2,1,pv(10));        % peak width
0485             qqu.ll=m(3,1,pv(10));       % peak slope on low side
0486             qqu.lh=m(4,1,pv(10));       % peak slope on high side
0487             qqu.gu=pact(2);             % guess rate
0488             qqu.la=pact(1);             % lapse rate
0489             sw=qqu.w*qqu.ll*qqu.lh/(qqu.ll+qqu.lh);   % normalized distance of inflections from peak
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,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)        % 'E/e': plot evolving Psychometric function on figure 103
0510         figure(pv(12)+3);
0511         psyevo=zeros(ntrial+1,nxevo);   % space for evolving pdf
0512         if unimode                      % unimodal psychometric function
0513             qqu.pk=m(1,jj,pv(10));      % peak position
0514             qqu.w=m(2,1,pv(10));        % peak width
0515             qqu.ll=m(3,1,pv(10));       % peak slope on low side
0516             qqu.lh=m(4,1,pv(10));       % peak slope on high side
0517             qqu.gu=pact(2);             % guess rate
0518             qqu.la=pact(1);             % lapse rate
0519             sw=qqu.w*qqu.ll*qqu.lh/(qqu.ll+qqu.lh);   % normalized distance of inflections from peak
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);     % peak position
0523                 qqu.w=mnt(2,1,pv(10),iet);        % peak width
0524                 qqu.ll=mnt(3,1,pv(10),iet);     % peak slope on low side
0525                 qqu.lh=mnt(4,1,pv(10),iet);     % peak slope on high side
0526                 sw=qqu.w*qqu.ll*qqu.lh/(qqu.ll+qqu.lh);   % normalized distance of inflections from peak
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'  % put trial 0 at the top
0533         else                            % monotonic psychometric function (not implemented)
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)   % 'V/v': plot convergence on figure 102
0540         figure(pv(12)+2);
0541         if unimode                      % unimodal psychometric function
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                            % monotonic psychometric function
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   % main loop for each probe value
0585 if pf(29) || pf(30)
0586     fclose(ofid);
0587 end

Generated on Tue 19-Sep-2017 12:07:31 by m2html © 2003