# v_psycdigit

## PURPOSE

V_PSYCDIGIT measures psychometric function using TIDIGITS stimuli

## SYNOPSIS

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

## DESCRIPTION

    V_PSYCDIGIT measures psychometric function using TIDIGITS stimuli

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

     Inputs:
             proc    processing function handle (e.g. @v_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 v_psycest or v_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:

- v_activlev V_ACTIVLEV Measure active speech level as in ITU-T P.56 [LEV,AF,FSO]=(sp,FS,MODE)
- v_psycest V_PSYCEST estimate multiple psychometric functions
- v_psycestu V_PSYCESTU estimate unimodal psychometric function
- v_psychofunc V_PSYCHOFUNC Calculate psychometric functions: trial success probability versus SNR
- v_readsph V_READSPH Read a SPHERE/TIMIT format sound file [Y,FS,WRD,PHN,FFX]=(FILENAME,MODE,NMAX,NSKIP)
- v_readwav V_READWAV Read a .WAV format sound file [Y,FS,WMODE,FIDX]=(FILENAME,MODE,NMAX,NSKIP)
- v_writewav V_WRITEWAV Creates .WAV format sound files FIDX=(D,FS,FILENAME,MODE,NSKIP,MASK,MAD)

This function is called by:

## SOURCE CODE

    0001 function [m,v]=v_psycdigit(proc,r,mode,p,q,xp,noise,fn,dfile,ofile)
    0002 %V_PSYCDIGIT measures psychometric function using TIDIGITS stimuli
    0003 %
    0004 % Usage:
    0005 %         (1)[m,v]=v_psycdigit([],[],'GMPWrn10',[],[],[],[],[],dfile);
    0006 %                       % measure SRT using addditive white noise, repetitions allowed, data in .WAV format
    0007 %         (2)[m,v]=v_psycdigit(@v_specsub,[],'GMPn10',[],[],[],[],[],dfile);
    0008 %                       % compare spectral subtraction with unprocessed noisy speech
    0009 %
    0010 % Inputs:
    0011 %         proc    processing function handle (e.g. @v_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 v_psycest or v_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: v_psycdigit.m 10865 2018-09-21 17:22:45Z 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 Lesser General Public License as published by
    0090 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0097 %
    0098 %   You can obtain a copy of the GNU Lesser General Public License from
    0099 %   https://www.gnu.org/licenses/ .
    0100 %
    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 v_psycest or v_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]=v_readwav([dfile tok{1,1}]); % get the first speech token to set fs
    0261 else
    0262     [s,fs]=v_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]=v_psycestu(-zmodel,p,q,xp); % initialize all models
    0277     [pact,qact]=v_psycestu(0);  % save the actual parameters
    0278 else
    0279     [xx,ii,m,v]=v_psycest(-zmodel,p,q,xp); % initialize all models
    0280     [pact,qact]=v_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]=v_readwav([dfile tok{isp,1}]); % get the speech token
    0357     else
    0358         [s,fs]=v_readsph([dfile tok{isp,1}]); % get the speech token
    0359     end
    0360     s=[zeros(pv(6)*round(fs/10),1); v_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                     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;   % 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]=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)       % '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             v_psycestu(jj);
    0467         else
    0468             v_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                 v_psycestu(jj);
    0475             else
    0476                 v_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,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)        % '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