function [ varargout ] = baumgartner2014_pmv2ppp( varargin )
%BAUMGARTNER2014_PMV2PPP Performance predictions from PMVs of baumgartner2014
% Usage: [ qe,pe,eb ] = baumgartner2014_pmv2ppp( p,tang,rang );
% [ qe,pe,eb ] = baumgartner2014_pmv2ppp( p,tang,rang,exptang );
% [ qe,pe,eb ] = baumgartner2014_pmv2ppp( pred );
%
% Input parameters:
% p : prediction matrix (response PMVs)
% tang : possible polar target angles. As default, ARI's MSP
% polar angles in the median SP is used.
% rang : polar angles of possible response angles.
% As default regular 5 deg.-sampling is used (-90:5:265).
% pred : structure including all input parameters above.
%
% Output parameters:
% qe : quadrant error rate
% pe : local polar RMS error in degrees
% eb : elevation bias in degrees; QEs and up-rear quadrant excluded
%
% BAUMGARTNER2014_PMV2PPP(...) retrieves commonly used PPPs (Psychoacoustic
% performance parameters) for sagittal-plane (SP) localization like quadrant
% error (QE), local polar RMS error (PE), and elevation bias (EB) from
% response PMVs (probability mass vectors) predicted by a localization
% model. PPPs are retreived either for a specific polar target angle or as
% an average across all available target angles. The latter is the
% default.
%
% BAUMGARTNER2014_PMV2PPP needs the following optional parameter in order
% to retrieve the PPPs for a specific (set of) target angles:
%
% 'exptang',exptang experimental polar target angles
%
% BAUMGARTNER2014_PMV2PPP accepts the following flag:
%
% 'print' Display the outcomes.
% 'QE_PE_EB' Compute QE, PE, and EB. This is the default.
% 'QE' Compute QE.
% 'PE' Compute PE.
% 'EB' Compute EB.
% 'absPE' Compute absolute polar error.
% 'chance' Compute chance performance for QE and PE.
%
% Example:
% ---------
%
% To evaluate chance performance of QE and PE use :
%
% [qe,pe] = baumgartner2014_pmv2ppp('chance');
%
% References:
% R. Baumgartner, P. Majdak, and B. Laback. Modeling sound-source
% localization in sagittal planes for human listeners. The Journal of the
% Acoustical Society of America, 136(2):791--802, 2014.
%
%
% Url: http://amtoolbox.org/amt-1.4.0/doc/modelstages/baumgartner2014_pmv2ppp.php
% #StatusDoc: Perfect
% #StatusCode: Perfect
% #Verification: Verified
% #Requirements: SOFA CircStat M-SIGNAL M-Stats O-Statistics
% #Author: Robert Baumgartner (2014), Acoustics Research Institute, Vienna, Austria
% This file is licensed unter the GNU General Public License (GPL) either
% version 3 of the license, or any later version as published by the Free Software
% Foundation. Details of the GPLv3 can be found in the AMT directory "licences" and
% at <https://www.gnu.org/licenses/gpl-3.0.html>.
% You can redistribute this file and/or modify it under the terms of the GPLv3.
% This file is distributed without any warranty; without even the implied warranty
% of merchantability or fitness for a particular purpose.
definput.import={'baumgartner2014_pmv2ppp'};
[flags,kv]=ltfatarghelper({'p','tang','rang','exptang'},definput,varargin);
if isstruct(kv.p) % input as *pred* structure
p = kv.p.p;
kv.rang = kv.p.rang;
kv.tang = kv.p.tang;
else
p = kv.p;
end
if flags.do_chance
p = ones(length(kv.rang),length(kv.tang));
end
if size(p,1) == 49 % rang: default for baumgartner2013
kv.rang=-30:5:210;
end
p = p./repmat(sum(p),length(kv.rang),1); % ensure probability mass vectors
tang = kv.tang(:);
rang = kv.rang(:);
nt = length(tang);
if not(flags.do_absPE)
qet = zeros(nt,1); % QE for each target angle
pet = zeros(nt,1); % PE for each target angle
ebt = zeros(nt,1); % EB for each target angle
isnotuprear = false(nt,1);
for ii = 1:nt % for all target positions
d = tang(ii)-rang; % wraped angular distance between tang & rang
iduw = (d < -180) | (180 < d); % 180deg-unwrap indices
d(iduw) = mod(d(iduw) + 180,360) - 180; % 180 deg unwrap
d = abs(d); % absolut distance
qet(ii) = sum( p(d>=90,ii) );
pc = p(d<90,ii); % pmv for conditional probability excluding QEs
pc = pc/sum(pc); % normalization to sum=1
pet(ii) = sqrt( sum( pc .* (d(d<90)).^2 )); % RMS of expected difference
if tang(ii) < 80
ebt(ii) = sum( pc .* rang(d<90) ) - tang(ii); % expectancy value of rang - tang
isnotuprear(ii) = true;
elseif tang(ii) > 180 % elevation instead of polar angle
ebt(ii) = -( sum( pc .* rang(d<90) ) - tang(ii) );
else % exclude up-rear quadrant
isnotuprear(ii) = false;
end
end
ebt = ebt(isnotuprear);
if ~isempty(kv.exptang)
qetb = (qet(1)+qet(end))/2; % boundaries for extang
petb = (pet(1)+pet(end))/2;
ebtb = (ebt(1)+ebt(end))/2;
extang = tang(:); % extended tang for targets outside
exqet = qet(:);
expet = pet(:);
expb = ebt(:);
if min(extang)>-90;
extang = [-90; extang];
exqet = [qetb; exqet];
expet = [petb; expet];
expb = [ebtb; expb];
isnotuprear = [true;isnotuprear];
end
if max(extang)<270;
extang = [extang; 270];
exqet = [exqet; qetb];
expet = [expet; petb];
expb = [expb; ebtb];
isnotuprear = [isnotuprear;true];
end
qet = interp1(extang,exqet,kv.exptang);
pet = interp1(extang,expet,kv.exptang);
excluderu = kv.exptang < 80 | kv.exptang > 180;
ebt = interp1(extang(isnotuprear),expb,kv.exptang(excluderu));
end
qe = mean(qet)*100;
pe = mean(pet);
eb = mean(ebt);
if isempty(flags.ppp) || flags.do_QE_PE_EB
varargout{1} = qe;
varargout{2} = pe;
varargout{3} = eb;
elseif flags.do_PE
varargout{1} = pe;
elseif flags.do_QE
varargout{1} = qe;
elseif flags.do_EB
varargout{1} = eb;
end
else % flags.do_absPE
apet = zeros(nt,1);
for ii = 1:nt % for all target positions
d = tang(ii)-rang; % wraped angular distance between tang & rang
iduw = (d < -180) | (180 < d); % 180-deg unwrap indices
d(iduw) = mod(d(iduw) + 180,360) - 180; % 180-deg unwrap
d = abs(d); % absolut distance
apet(ii) = sum( p(:,ii) .* d); % absolute polar angle error for target ii
end
ape = mean(apet);
varargout{1} = ape;
end
%% Print Output
if flags.do_print
if isempty(flags.ppp) || flags.do_QE_PE_EB
amt_disp(sprintf('Quadrant errors (%%) \t\t %4.1f',qe),'documentation');
amt_disp(sprintf('Local polar RMS error (deg) \t %4.1f',pe),'documentation');
if nargout==3
amt_disp(sprintf('Local polar bias (deg) \t\t %4.1f',eb),'documentation');
end
elseif flags.do_PE
amt_disp(sprintf('Quadrant errors (%%) \t\t %4.1f',qe),'documentation');
elseif flags.do_QE
amt_disp(sprintf('Local polar RMS error (deg) \t %4.1f',pe),'documentation');
elseif flags.do_EB
amt_disp(sprintf('Local elevation bias (deg) \t\t %4.1f',eb),'documentation');
elseif flags.do_absPE
amt_disp(sprintf('Absolute polar error (deg) \t\t %4.1f',ape),'documentation');
end
end