0001 function [xx,ii,m,v]=v_psycestu(iq,x,r,xp)
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 persistent pp qq fl fh fxs tx ts tl th
0038 if iq<0
0039 if iq~=-1
0040 error('Cannot yet have multiple models');
0041 end
0042 pp=x;
0043 qq=r;
0044 pp(7:2:9)=max(pp(7:2:9),qq.sl);
0045 end
0046 nxslh=qq.nxslh;
0047 nx=nxslh(1);
0048 ns=nxslh(2);
0049 nl=nxslh(3);
0050 nh=nxslh(4);
0051 la=pp(1);
0052 gu=pp(2);
0053 if iq<0
0054
0055 fl=ones(nx,nl,ns);
0056 fh=ones(nx,nh,ns);
0057 fxs=ones(nx,1,ns);
0058
0059 tx=linspace(pp(3),pp(4),nx)';
0060 ts=linspace(pp(5),pp(6),ns)';
0061 tl=linspace(log(pp(7)),log(pp(8)),nl)';
0062 th=linspace(log(pp(9)),log(pp(10)),nh)';
0063 elseif iq>0 && nargin==3
0064 if iq~=1
0065 error('Cannot yet have multiple models');
0066 end
0067 r0=r==0;
0068
0069 mskx=tx>x;
0070 nxm=sum(mskx);
0071 ets=reshape(exp(-ts),[1 1 ns]);
0072 sd0=1+ets.^2;
0073 sm0=2*ets;
0074 if any(mskx)
0075 fl(mskx,:,:)=fl(mskx,:,:).*(r0+(1-2*r0)*(gu+(1-gu-la)*(repmat(sd0,[nxm,nl,1])+repmat(sm0,[nxm,nl,1]).*repmat(cosh((tx(mskx)-x)*exp(tl')),[1 1 ns])).^(-1)));
0076 end
0077 if ~all(mskx)
0078 fh(~mskx,:,:)=fh(~mskx,:,:).*(r0+(1-2*r0)*(gu+(1-gu-la)*(repmat(sd0,[nx-nxm,nh,1])+repmat(sm0,[nx-nxm,nh,1]).*repmat(cosh((x-tx(~mskx))*exp(th')),[1 1 ns])).^(-1)));
0079 end
0080 else
0081 xx=pp;
0082 ii=qq;
0083 end
0084
0085
0086
0087 sfl=sum(fl,2);
0088 fl=fl./sfl(:,ones(nl,1),:);
0089 sfh=sum(fh,2);
0090 fh=fh./sfh(:,ones(nh,1),:);
0091 fxs=fxs.*sfl.*sfh;
0092 fxs=fxs/sum(fxs(:));
0093
0094
0095
0096 px=squeeze(sum(fxs,3));
0097 ps=squeeze(sum(fxs,1));
0098 pl=reshape(permute(fl,[2 1 3]),nl,nx*ns)*fxs(:);
0099 ph=reshape(permute(fh,[2 1 3]),nh,nx*ns)*fxs(:);
0100
0101
0102
0103 m=zeros(4,1,3);
0104 m(:,1,1)=[tx'*px ts'*ps tl'*pl th'*ph]';
0105 vraw=[tx'.^2*px ts'.^2*ps tl'.^2*pl th'.^2*ph]'-m(:,1,1).^2;
0106
0107
0108 [flm,iflm]=max(fl,[],2);
0109 [fhm,ifhm]=max(fh,[],2);
0110 fxsm=fxs.*flm.*fhm;
0111 [pxpk,i]=max(fxsm(:));
0112 j=1+floor((i-1)/nx);
0113 i=i-nx*(j-1);
0114 jl=iflm(i,1,j);
0115 jh=ifhm(i,1,j);
0116
0117
0118
0119
0120 m(:,1,2)=[tx(i) ts(j) tl(jl) th(jh)]';
0121
0122 [pxpk,i]=max(px);
0123 if i>1 && i<nx
0124 [dum,j]=v_quadpeak(px(i-1:i+1));
0125 i=i+j-2;
0126 end
0127 xmm=(2-i)*tx(1)+(i-1)*tx(2);
0128
0129 [pxpk,i]=max(ps);
0130 if i>1 && i<ns
0131 [dum,j]=v_quadpeak(ps(i-1:i+1));
0132 i=i+j-2;
0133 end
0134 smm=(2-i)*ts(1)+(i-1)*ts(2);
0135
0136 [pxpk,i]=max(pl);
0137 if i>1 && i<nl
0138 [dum,j]=v_quadpeak(pl(i-1:i+1));
0139 i=i+j-2;
0140 end
0141 lmm=(2-i)*tl(1)+(i-1)*tl(2);
0142
0143 [pxpk,i]=max(ph);
0144 if i>1 && i<nh
0145 [dum,j]=v_quadpeak(ph(i-1:i+1));
0146 i=i+j-2;
0147 end
0148 hmm=(2-i)*th(1)+(i-1)*th(2);
0149 m(:,1,3)=[xmm smm lmm hmm]';
0150 m(3:4,:,:)=exp(m(3:4,:,:));
0151 mfact=(m(3,1,:)+m(4,1,:))./(m(3,1,:).*m(4,1,:));
0152 m(2,1,:)=m(2,1,:).*mfact;
0153 v=[vraw(1); vraw(2)*mfact(1).^2; vraw(3:4).*m(3:4,1,1).^2];
0154
0155
0156
0157
0158
0159 if ~nargout && iq==1
0160 subplot(121)
0161 imagesc(tx,ts,squeeze(fxs)');
0162 axis 'xy'
0163 xlabel('Peak position (dB)');
0164 ylabel('Normalized semi-width (dB)');
0165 subplot(122)
0166 imagesc(th,tl,reshape(permute(fl,[2 1 3]),nl,nx*ns)*reshape(permute(fh.*fxs(:,ones(nh,1),:),[1 3 2]),nx*ns,nh));
0167 axis 'xy'
0168 xlabel('Ln down slope (ln prob/dB)');
0169 ylabel('Ln up slope (ln prob/dB)');
0170 end
0171
0172
0173 if iq~=0
0174 xt=tx;
0175 nxt=numel(xt);
0176 ets=reshape(exp(-ts),[1 1 ns]);
0177 sd0=1+ets.^2;
0178 sm0=2*ets;
0179 hh=zeros(nxt,1);
0180 for i=1:nxt
0181 y=xt(i);
0182 mskx=tx>y;
0183 nxm=sum(mskx);
0184 flm=gu+(1-gu-la)*(repmat(sd0,[nxm,nl,1])+repmat(sm0,[nxm,nl,1]).*repmat(cosh((tx(mskx)-y)*exp(tl')),[1 1 ns])).^(-1);
0185 fhm=gu+(1-gu-la)*(repmat(sd0,[nx-nxm,nh,1])+repmat(sm0,[nx-nxm,nh,1]).*repmat(cosh((y-tx(~mskx))*exp(th')),[1 1 ns])).^(-1);
0186
0187 fl1=fl;
0188 fh1=fh;
0189 fl1(mskx,:,:)=fl1(mskx,:,:).*flm;
0190 fh1(~mskx,:,:)=fh1(~mskx,:,:).*fhm;
0191 sfl1=sum(fl1,2);
0192 fl1=fl1./sfl1(:,ones(nl,1),:);
0193 sfh1=sum(fh1,2);
0194 fh1=fh1./sfh1(:,ones(nh,1),:);
0195 fxs1=fxs.*sfl1.*sfh1;
0196 pr=sum(fxs1(:));
0197 fxs1=fxs1/pr;
0198 px1=squeeze(sum(fxs1,3));
0199 ps1=squeeze(sum(fxs1,1));
0200 pl1=reshape(permute(fl1,[2 1 3]),nl,nx*ns)*fxs1(:);
0201 ph1=reshape(permute(fh1,[2 1 3]),nh,nx*ns)*fxs1(:);
0202
0203 fl0=fl;
0204 fh0=fh;
0205 fl0(mskx,:,:)=fl0(mskx,:,:).*(1-flm);
0206 fh0(~mskx,:,:)=fh0(~mskx,:,:).*(1-fhm);
0207 sfl0=sum(fl0,2);
0208 fl0=fl0./sfl0(:,ones(nl,1),:);
0209 sfh0=sum(fh0,2);
0210 fh0=fh0./sfh0(:,ones(nh,1),:);
0211 fxs0=fxs.*sfl0.*sfh0;
0212 fxs0=fxs0/(1-pr);
0213 px0=squeeze(sum(fxs0,3));
0214 ps0=squeeze(sum(fxs0,1));
0215 pl0=reshape(permute(fl0,[2 1 3]),nl,nx*ns)*fxs0(:);
0216 ph0=reshape(permute(fh0,[2 1 3]),nh,nx*ns)*fxs0(:);
0217
0218 h1=[(tx(1)-tx(2))*px1'*log(px1) (ts(1)-ts(2))*ps1'*log(ps1) (tl(1)-tl(2))*pl1'*log(pl1) (th(1)-th(2))*ph1'*log(ph1)]';
0219 h0=[(tx(1)-tx(2))*px0'*log(px0) (ts(1)-ts(2))*ps0'*log(ps0) (tl(1)-tl(2))*pl0'*log(pl0) (th(1)-th(2))*ph0'*log(ph0)]';
0220 hh(i)=qq.cs*(pr*h1+(1-pr)*h0);
0221 end
0222 [hmin,ih]=min(hh);
0223 xx=xt(ih);
0224 ii=1;
0225 end
0226
0227
0228
0229
0230
0231
0232
0233
0234
0235 plow=max(min(0.1*gu,0.001),1e-6);
0236
0237 pcx=cumsum(px);
0238 tmp=linspace(tx(max(find(pcx>plow,1)-1,1)),tx(find(pcx>(1-plow),1)),nx)';
0239 fxs=linres(fxs,1,tx,tmp);
0240 fl=linres(fl,1,tx,tmp);
0241 fh=linres(fh,1,tx,tmp);
0242 tx=tmp;
0243
0244 pcx=cumsum(ps);
0245 tmp=linspace(ts(max(find(pcx>plow,1)-1,1)),ts(find(pcx>(1-plow),1)),ns)';
0246 fxs=linres(fxs,3,ts,tmp);
0247 fl=linres(fl,3,ts,tmp);
0248 fh=linres(fh,3,ts,tmp);
0249 ts=tmp;
0250
0251
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261 function y=linres(x,d,a,b)
0262
0263 sz=size(x);
0264 n=sz(d);
0265 p=[d:numel(sz) 1:d-1];
0266 sz2=prod(sz)/n;
0267 z=reshape(permute(x,p),n,sz2);
0268 na=numel(a);
0269 nb=numel(b);
0270 [xx,ix]=sort([a(:); b(:)]);
0271 jx=zeros(1,na+nb);
0272 jx(ix)=(1:na+nb);
0273
0274 jb=jx(na+1:end)-(1:nb);
0275 y=zeros(nb,size(z,2));
0276 y(jb==0,:)=z(ones(sum(jb==0),1),:);
0277 y(jb==na,:)=z(na*ones(sum(jb==na),1),:);
0278 mk=(jb>0) & (jb<na);
0279 jmk=jb(mk);
0280 f=((b(mk)-a(jmk))./(a(1+jmk)-a(jmk)));
0281 fc=1-f;
0282 y(mk,:)=z(jmk,:).*fc(:,ones(1,sz2))+z(1+jmk,:).*f(:,ones(1,sz2));
0283 sz(d)=nb;
0284 y=ipermute(reshape(y,sz(p)),p);
0285
0286
0287
0288
0289
0290