THE AUDITORY MODELING TOOLBOX

Applies to version: 1.6.0

View the help

Go to function

localizationerror
Compute psychoacoustic performance parameters from a sound localization experiment

Program code:

function [varargout] = localizationerror(m,varargin)
%localizationerror Compute psychoacoustic performance parameters from a sound localization experiment
%   Usage: [accL, precL, accP, precP, querr] = localizationerror(m)
%          [res, meta, par] = localizationerror(m,errorflag)
%          [res, meta] = localizationerror(m,f,r,'perMacpherson2003');
%
%   Input parameters:
%     m  : Item list from a localization experiment where each row
%          is response from a trial and the columns are as follows:
%
%          - 1: Azimuth angle of the target (in degrees).
%
%          - 2: Elevation angle of the target (in degrees).
%
%          - 3: Azimuth angle of the response (in degrees).
%
%          - 4: Elevation angle of the response (in degrees).
%
%          - 5: Lateral angle of the target (in degrees).
%
%          - 6: Polar angle of the target (in degrees).
%
%          - 7: Lateral angle of the response (in degrees).
%
%          - 8: Polar angle of the response (in degrees).
%
%          - 9 to 11: Cartesian X, Y, and Z coordinates of the target. Required only for errorflags being
%            'rFBC' or 'rF2B'.
%
%          - 12 to 14: Cartesian X, Y, and Z coordinates of the response. Required only for errorflags being
%            'rFBC' or 'rF2B'.
%
%     errorflag : Error metric to be calculated.
%
%     f  : Regression statistics for the frontal targets obtained from [f,r] = LOCALIZATIONERROR(m, 'sirp');.
%          Required only for the errorflag being 'perMacpherson2003'.
%
%     r  : Regression statistics for the rear targets obtained from [f,r] = LOCALIZATIONERROR(m, 'sirp');.
%          Required only for the errorflag being 'perMacpherson2003'.
%
%   Output parameters:
%     accL  : Accuracy bias (in degrees) in the lateral dimension.
%     precL : Precision error (in degrees) in the lateral dimension.
%     accP  : Accuracy bias (in degrees) in the polar dimension.
%     precP : Precision error (in degrees) in the polar dimension.
%     querr : Quadrant error rate (in percentage).
%     res   : Error metric as requested by errorflag.
%     meta  : Additional metadata as requested by errorflag.
%     par   : Additional calculation as requested by errorflag.
%
%   LOCALIZATIONERROR(..) calculates various psychoacoustic performance
%   parameters from responses of a sound-localization experiment.
%
%   [accL, precL, accP, precP, querr] = LOCALIZATIONERROR(m) returns the metrics 'accL', 'precL',
%   'precP', and 'querr'.
%
%   [res, meta, par] = LOCALIZATIONERROR(m, errorflag) returns in res the error
%   metric requested by errorflag, with additional metadata and calculation results
%   provided in meta and par, respectively.
%
%   The errorflag can be one of the following:
%
%     'absaccE'            Absolute elevation accuracy bias (in degrees), i.e., the absolute value
%                          of the circular average of elevation errors. See magnitude of elevation
%                          bias in Middlebrooks (1999).
%
%     'absaccL'            Absolute lateral accuracy bias (in degrees), i.e., the absolute value
%                          of the average of the (signed) errors in the lateral dimension. See magnitude
%                          of lateral bias in Middlebrooks (1999).
%
%     'accabsL'            Averaged absolute lateral erorr (in degrees), i.e., the
%                          average of the absolute (unsigned) lateral errors. Reference unclear. Note that
%                          the naming as accuracy is confusing.
%
%     'accabsP'            Averaged absolute polar error (in degrees), i.e., the circular
%                          average of the absolute (unsigned) polar errors. Reference unclear. Note that the
%                          naming as accuracy is confusing.
%
%     'accE'               Elevation accuracy bias (in degrees), i.e., the circular average of
%                          elevation errors. Reference unclear.
%
%     'accL'               Lateral acuracy bias (in degrees), i.e., the average of lateral errors. See lateral
%                          bias in Majdak et al. (2010).
%
%     'accP'               Polar accuracy bias (in degrees), i.e., the circular average of polar errors.
%                          Reference unclear.
%
%     'accPnoquerr'        As 'accP' but with quadrant errors (see 'querr') removed. See local polar bias*
%                          in Majdak et al. (2010).
%
%     'corrcoefL'          Lateral correlation coeffficient, i.e., the correlation coefficient between the
%                          response and target angles in the lateral dimension. meta provides the p-value
%                          of the correlation significance. See lateral target-response correlation coefficient*
%                          in Majdak et al. (2011).
%
%     'corrcoefP'          Polar median-plane correlation coeffficient, i.e., the correlation coefficient between
%                          the response and target angles in the polar dimension, considering responses around the
%                          median plane only, i.e., within the lateral angles of +-30 ^circ. meta provides
%                          the p-value of the correlation significance. See polar target-response correlation coefficient*
%                          in Majdak et al. (2011).
%
%     'gainL'              Lateral gain. It uses 'gainLstats'. Reference unclear.
%
%     'gainLstats'         Statistic details from the calculation of the lateral gain. See regress for a detailed
%                          description of the structure. Reference unclear.
%
%     'gainP'              Polar gain (between -1 and 1) as an average of 'gainPfront' and 'gainPrear'. See
%                          polar gain in Macpherson and Middlebrooks (2000).
%
%     'gainPfront'         Frontal polar gain (between -1 and 1), i.e., the  of SIRP done for frontal targets. It uses
%                          'sirp'. See polar gain in Macpherson and Middlebrooks (2000).
%
%     'gainPrear'          Rear polar gain (between -1 and 1), i.e., the slope of SIRP done for rear targets. It uses
%                          'sirp'. See polar gain in Macpherson and Middlebrooks (2000).
%
%     'perMacpherson2003'  Polar error rate (in percentage) as the rate of deviations larger than 45^circ from the
%                          linear regression calculated by SIRP for polar angles. It requires the results from 'sirp'
%                          provided as a separate inputs (see Usage above).
%
%     'precE'              Elevation precision error (in degrees) i.e., circular standard deviation of elevation errors.
%                          Reference unclear.
%
%     'precL'              Lateral precision error (in degrees), i.e., circular standard deviation of lateral errors.
%                          See lateral precision error in Majdak et al. (2010).
%
%     'precLcentral'       As 'precL' but considering only responses +-30 ^circ
%                          in elevation around the horizontal meridian plane. Reference unclear.
%
%     'precLregress'       Lateral scatter (in degrees) around the lateral regression line,
%                          i.e., the RMS of the response deviation from the lateral regression
%                          line, considering only the quasi-veridical responses, i.e., deviations
%                          smaller then 45^circ. It uses 'gainLstats'. See lateral angle scatter*
%                          in Macpherson and Middlebrooks (2000).
%                          Note that the naming as precision is confusing.
%
%     'precP'              Polar precision error (in degrees), i.e., circular standard deviation of polar errors.
%                          Reference unclear.
%
%     'precPmedian'        As 'precP' but considering targets around the median plane, i.e., with lateral angles
%                          in the range of +-30 ^circ, only. Reference unclear.
%
%     'precPmedianlocal'   As 'precpPmedian' but  with quadrant errors (see 'querrMiddlebrooks') removed.
%                          Similar but not the same as rms local polar error in Middlebrooks (1999).
%
%     'precPnoquerr'       As 'precP' but with quadrant errors (see 'querr') removed. See local polar
%                          precision error in Majdak et al. (2010).
%
%     'precPregressFront'  Polar frontal scatter (in degrees) around the polar regression line,
%                          i.e., the RMS of the response deviation from the polar regression
%                          line calculated by SIRP considering only targets located in the front
%                          and only the quasi-veridical responses, i.e., deviations
%                          smaller then 45^circ and. It uses 'sirp'. See polar angle scatter*
%                          in Macpherson and Middlebrooks (2000).
%                          Note that the naming as precision is confusing.
%
%     'precPregressRear'   Polar rear scatter (in degrees). As 'precPregressFront' but calculated
%                          considering only targets located in the rear. It uses 'sirp'.
%                          See polar angle scatter in Macpherson and Middlebrooks (2000).
%                          Note that the naming as precision is confusing.
%
%     'precPregress'       Polar scatter (in degrees) as an average of 'precPregressFront' and
%                          'precPregressRear'. See polar angle scatter in Macpherson and Middlebrooks (2000).
%                          Note that the naming as precision is confusing.
%
%     'pVeridicalL'        Proportion (in percentage) of lateral quasi-verdical responses. It uses 'gainLstats'.
%                          Reference unclear.
%
%     'pVeridicalPfront'   Proportion (in percentage) of polar frontal quasi-verdical polar responses. See
%                          quasi-veridical responses in Macpherson and Middlebrooks (2000).
%
%     'pVeridicalPrear'    Proportion (in percentage) of polar rear quasi-verdical polar responses. See
%                          quasi-veridical responses in Macpherson and Middlebrooks (2000).
%
%     'pVeridicalP'        Proportion (in percentage) of polar quasi-verdical polar responses.See quasi-veridical
%                          responses in Macpherson and Middlebrooks (2000).
%
%     'querr'              Quadrant error rate (in percentage), i.e., the rate of responses with the weighted
%                          polar errors outside the quadrant defined by errors within +-45 ^circ. The
%                          weighting is done as a function the lateral angle of the target. meta and par*
%                          provide the information about the number of responses in the correct quadrants, the
%                          number of responses within the lateral range, and the chance rate to achieve
%                          the quadrant errors by chance. See quadrant errors in Majdak et al. (2010).
%
%     'querrMiddlebrooks'  Quadrant error rate (in percentage), i.e., the rate of responses with
%                          the unweighted polar errors outside the quadrant defined by absolute
%                          polar errors within +-90 ^circ. meta and par provide the information
%                          about the number of confusions and the number of responses within the lateral range.
%                          See quadrant errors in Middlebrooks (1999).
%
%     'rmsL'               Lateral RMS error (in degrees), i.e., the root of the average of squared
%                          lateral errors. See rms lateral error in Middlebrooks (1999).
%
%     'rmsPmedian'         Polar median-plane RMS error (in degrees), i.e., the root of the average of squared
%                          polar errors considering only targets around the median plane, i.e., with lateral
%                          angles in the range of +-30 ^circ, only.  Reference unclear.
%
%     'rmsPmedianlocal'    As 'rmsPmedian' but excluding responses with quadrant errors as in
%                          'querrMiddlebrooks', i.e., absolute polar errors larger than 90 ^circ.
%                          See rms local polar error in Middlebrooks (1999).
%
%     'SCC'                Spherical correlation coefficient. See SCC in Carlile et al. (1997).
%
%     'sirp'               Polar regression statistics by means of the selective iterative regression procedure
%                          (SIRP) to exclude outliers and reversals from the computation of the regression lines
%                          used to calculate gain-based metrics. The SIRP is run 100 times and the result
%                          including the most responses (in other words, excluding the least outliers) is selected.
%                          Only targets around the median plane are considered, i.e., within the lateral range of
%                          +-30 ^circ. Outlier distance criterion is 40^circ. The SIRP is run separately
%                          for targets located in the frontal and the rear hemisphere with the results stored in res*
%                          and meta respectively (see regress for a detailed description of the structure).
%                          See p. 1873 in Macpherson and Middlebrooks (2000).
%
%     'rFBC'               Rate of front-back confusion (in percentage) in terms of being in the incorrect hemifield
%                          separated by the frontal plane. Responses on the interaural axis are not considered.
%                          Requires additional columns in the input. See front-back confusions in Carlile et al.
%                          (1997). Requires additional columns in the input.
%
%     'rF2B'               As 'rFBC' but considering front-to-back confusions only. Requires additional columns
%                          in the input. See McLachlan et al. (2024).
%
%     'rUDC'               Rate of up-down confusion (in percentage) in terms of being in the incorrect hemifield
%                          separated by the horizontal meridian plane. Responses directly on that plane are not
%                          considered. Requires additional columns in the input. See up-down confusions in
%                          Carlile et al. (1997) and McLachlan et al. (2024).
%
%     'slopePfront'        Polar frontal slope angle (in degrees) calcalated from the slope of the regression
%                          line for the frontal targets. Uses 'gainPfront'. Reference unclear.
%
%     'slopePrear'         Polar rear slope angle (in degrees) calcalated from the slope of the regression line
%                          for the rear targets. Uses 'gainPrear'. Reference unclear.
%
%     'slopeP'             Polar slope angle (in degrees) calcalated from the slope of the regression line for
%                          all targets. Uses 'gainP'. Reference unclear.
%
%
%
%
%
%   References:
%     G. McLachlan, P. Majdak, J. Reijniers, M. Mihocic, and H. Peremans.
%     Insights into dynamic sound localisation: A direction-dependent
%     comparison between human listeners and a Bayesian model, Apr. 2024.
%     
%     E. A. Macpherson and J. C. Middlebrooks. Localization of brief sounds:
%     Effects of level and background noise. The Journal of the Acoustical
%     Society of America, 108(1834), 2000.
%     
%     P. Majdak, M. J. Goupell, and B. Laback. Two-dimensional localization
%     of virtual sound sources in cochlear-implant listeners. Ear Hear,
%     32:198--208, 2011.
%     
%     P. Majdak, M. J. Goupell, and B. Laback. 3-D localization of virtual
%     sound sources: Effects of visual environment, pointing method and
%     training. Atten Percept Psycho, 72:454--469, 2010.
%     
%     J. C. Middlebrooks. Virtual localization improved by scaling
%     nonindividualized external-ear transfer functions in frequency. The
%     Journal of the Acoustical Society of America, 106:1493--1510, 1999.
%     
%     E. A. Macpherson and J. C. Middlebrooks. Vertical-plane sound
%     localization probed with ripple-spectrum noise. The Journal of the
%     Acoustical Society of America, 114:430--445, 2003.
%     
%     S. Carlile, P. Leong, and S. Hyams. The nature and distribution of
%     errors in sound localization by human listeners. Hearing research,
%     114(1):179--196, 1997.
%     
%
%   See also: exp_baumgartner2014 exp_barumerli2023
%
%   Url: http://amtoolbox.org/amt-1.6.0/doc/common/localizationerror.php


%   #Author: Piotr Majdak (2021): original implementation
%   #Author: Robert Baumgartner (2021): various additions, especially those based on sirpMacpherson2000
%   #Author: Glen McLachlan (2023): addition of rFBC, rF2B, rUDC, and maePnorev
%   #Author: Piotr Majdak (2024): major rework of the documentation for AMT 1.6, sirpMacpherson2000 renamed to sirp. maePnorev removed.
%   #Author: Michael Mihocic (2024): 360-deg wrapping fixed to be compatible to Octave code


% 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.

%% Check input options

definput.import={'localizationerror'};

[flags,kv]=ltfatarghelper({'f','r'},definput,varargin);


if size(m,2)==1,   m=m'; end


if isempty(flags.errorflag)
    nargout=6;
    thr=45;

    accL=mean(m(:,7)-m(:,5));
    varargout{1}=accL;
    precL=sqrt(sum((m(:,7)-m(:,5)-accL).^2)/size(m,1));
    varargout{2}=precL;
    w=0.5*cos(deg2rad(m(:,7))*2)+0.5;
    varargout{3}=mynpi2pi(circmean(m(:,8)-m(:,6),w));
    [x,mag]=circmean(m(:,8)-m(:,6));
    precP=rad2deg(acos(mag)*2)/2;
    varargout{4}=precP;
    idx=find((abs(mynpi2pi(m(:,8)-m(:,6))).*w)<=thr);
    querr=100-size(idx,1)/size(m,1)*100;
    varargout{5}=querr;
else
    %   if iscell(errorflag), errorflag=errorflag{1}; end
    %   if ~ischar(errorflag)
    %     error('ERR must be a string');
    %   end
    %   if size(errorflag,1)~=1, errorflag=errorflag'; end
    nargout=3;
    metadata=[];
    par=[];
    switch flags.errorflag
        case 'accL'
            accL=mean(m(:,7)-m(:,5));
            varargout{1}=accL;
            meta.ylabel='Lateral bias (deg)';
        case 'absaccL'
            accL=mean(m(:,7)-m(:,5));
            varargout{1}=abs(accL);
            meta.ylabel='Lateral absolute bias (deg)';
        case 'accabsL' %i.e. mean absolute error
            accL=mean(abs(m(:,7)-m(:,5)));
            varargout{1}=accL;
            meta.ylabel='Lateral absolute error (deg)';
        case 'corrcoefL'
            [r,p] = corrcoef([m(:,5) m(:,7)]);
            varargout{1}=r(1,2);
            par.par1=p(1,2);
            meta.ylabel='Lateral correlation coefficient';
            meta.par1 = 'p for testing the hypothesis of no correlation';
        case 'corrcoefP'
            idx=find(abs(m(:,7))<=30);
            [r,p] = corrcoef([m(idx,6) m(idx,8)]);
            varargout{1}=r(1,2);
            par.par1=p(1,2);
            meta.ylabel='Polar correlation coefficient';
            meta.par1 = 'p for testing the hypothesis of no correlation';
        case 'precL'
            accL=mean(m(:,7)-m(:,5));
            precL=sqrt(sum((m(:,7)-m(:,5)-accL).^2)/size(m,1));
            varargout{1}=precL;
            meta.ylabel='Lateral precision error (deg)';
        case 'precLcentral'
            idx=find(abs(m(:,4))<=30); % responses close to horizontal meridian
            accL=mean(m(idx,7)-m(idx,5));
            precL=sqrt(sum((m(idx,7)-m(idx,5)-accL).^2)/length(idx));
            varargout{1}=precL;
            meta.ylabel='Central lateral precision error (deg)';
        case 'rmsL'
            idx=find(abs(m(:,7))<=60);
            m=m(idx,:);
            rmsL=sqrt(sum((mynpi2pi(m(:,7)-m(:,5))).^2)/size(m,1));
            varargout{1}=rmsL;
            meta.ylabel='Lateral RMS error (deg)';
        case 'accP'
            varargout{1}=mynpi2pi(circmean(m(:,8)-m(:,6)));
            meta.ylabel='Polar bias (deg)';
        case 'accabsP'
            varargout{1}=mynpi2pi(circmean(abs(m(:,8)-m(:,6))));
            meta.ylabel='Polar absolute bias (deg)';
        case 'accPnoquerr'
            w=0.5*cos(deg2rad(m(:,7))*2)+0.5;
            idx=find(w>0.5);
            m=m(idx,:); w=w(idx,:);
            idx=find((abs(mynpi2pi(m(:,8)-m(:,6))).*w)<=45);
            varargout{1}=mynpi2pi(circmean(m(idx,8)-m(idx,6)));
            par.par1=length(idx);
            meta.ylabel='Local polar bias (deg)';
            meta.par1='Number of considered responses';
        case 'precP'
            w=0.5*cos(deg2rad(m(:,7))*2)+0.5;
            precP=circstd(m(:,8)-m(:,6),w);
            varargout{1}=precP;
            meta.ylabel='Polar precision (deg)';
        case 'precPnoquerr'
            %         w=0.5*cos(deg2rad(m(:,7))*2)+0.5;
            %         idx=find(w>0.5);
            %         thr=45; % threshold for the choice: quadrant error or not
            %         m=m(idx,:); w=w(idx,:);
            range=360; % polar range of targets
            thr=45; % threshold for the choice: quadrant error or not
            latmax=0.5*180*(acos(4*thr/range-1))/pi;
            idx=find(abs(m(:,7))<latmax); % targets outside not considered
            m=m(idx,:);
            w=0.5*cos((m(:,7))*2*pi/180)+0.5;
            idx=find((abs(mynpi2pi(m(:,8)-m(:,6))).*w)<=thr);
            m=m(idx,:);
            w=0.5*cos((m(:,7))*pi/180*2)+0.5;
            precP=circstd(m(:,8)-m(:,6),w);
            varargout{1}=real(precP);
            meta.ylabel='Local polar precision error (deg)';
        case 'querr'
            range=360; % polar range of targets
            thr=45; % threshold for the choice: quadrant error or not
            latmax=0.5*180*(acos(4*thr/range-1))/pi;
            idx=find(abs(m(:,7))<latmax); % targets outside not considered
            m=m(idx,:);
            w=0.5*cos(pi*(m(:,7))/180*2)+0.5; % weigthing of the errors
            corrects=size(find((abs(mynpi2pi(m(:,8)-m(:,6))).*w)<=thr),1);
            totals=size(m,1);
            chancerate=2*thr/range./w;
            wrongrate=100-mean(chancerate)*100;
            querr=100-corrects/totals*100;
            varargout{1}=querr;
            meta.ylabel='Weighted quadrant errors (%)';
            par.par1=corrects;
            meta.par1='Number of correct responses';
            par.par2=totals;
            meta.par2='Number of responses within the lateral range';
            par.par3=wrongrate; % valid only for targets between 0 and 180°
            meta.par3='Wrong rate, valid only for targets between 0 and 180°';

        case 'querrMiddlebrooks' % as in Middlebrooks (1999); comparable to COC from Best et al. (2005)
            idx=abs(m(:,7))<=30;
            m=m(idx,:);
            idx=find((abs(mynpi2pi(m(:,8)-m(:,6))))>90);
            querr=size(idx,1)/size(m,1)*100;
            % chance performance (FIXME: not working with exp_baumgartner2014!)
            %       tang = m(:,6);
            %       rang = min(tang):max(tang); % assume listener knows target range
            %       p = ones(length(rang),length(tang)); % equally distributed response probabilities
            %       chance = baumgartner2014_pmv2ppp(p,tang,rang,'QE');
            varargout{1}=querr;
            meta.ylabel='Quadrant errors (%)';
            par.par1=size(idx,1);
            meta.par1='Number of confusions';
            par.par2=size(m,1);
            meta.par2='Number of responses within the lateral range';
            %       par.chance=chance;
            %       meta.chance='Chance performance';

        case 'querrMiddlebrooksRAU' % as in Middlebrooks (1999) but calculated in rau
            idx=find(abs(m(:,7))<=30);
            m=m(idx,:);
            idx=find((abs(mynpi2pi(m(:,8)-m(:,6))))<=90);
            querr=100-size(idx,1)/size(m,1)*100;
            varargout{1}=rau(querr,size(m,1));
            meta.ylabel='RAU quadrant errors';
            par.par1=size(idx,1);
            meta.par1='Number of confusions';
            par.par2=size(m,1);
            meta.par2='Number of responses within the lateral range';
        case 'rFBC'
            idx1 = find(sign(m(:,9))~=sign(m(:,12)) & abs(m(:,9))>0.01);
            idx2 = abs(m(:,9))>0.01;
            varargout{1}=size(idx1,1)/size(idx2,1)*100;
            meta.ylabel='FBC rate (%)';
        case 'rF2B' %front-to-back
            idx1 = find(sign(m(:,9))~=sign(m(:,12)) & abs(m(:,9))>0.01);
            idx2 = find(sign(m(:,9))~=sign(m(:,12)) & (m(:,9))>0.01);
            varargout{1}=size(idx2,1)/size(idx1,1)*100;
            meta.ylabel='F2B rate (%)';
        case 'rUDC'
            idx1 = find(sign(m(:,11))~=sign(m(:,14)) & abs(m(:,11))>0.01);
            idx2=abs(m(:,11))>0.01;
            varargout{1}=size(idx1,1)/size(idx2,1)*100;
            meta.ylabel='UDC rate (%)';
        case 'accE'
            varargout{1}=mynpi2pi(circmean(m(:,4)-m(:,2)));
            meta.ylabel='Elevation bias (deg)';
        case 'absaccE'
            varargout{1}=abs(mynpi2pi(circmean(m(:,4)-m(:,2))));
            meta.ylabel='Absolute elevation bias (deg)';
        case 'precE'
            precE=circstd(m(:,4)-m(:,2));
            varargout{1}=precE;
            meta.ylabel='Elevation precision error (deg)';
        case 'lateral'        % obsolete
            nargout=2;
            varargout{1}=accL;
            varargout{2}=precL;
        case 'polar'          % obsolete
            nargout=3;
            varargout{1}=accP;
            varargout{2}=precP;
            varargout{3}=querr;
        case 'precPmedianlocal'
            idx=find(abs(m(:,7))<=30);
            m=m(idx,:);
            idx=find(abs(mynpi2pi(m(:,8)-m(:,6)))<90);
            if isempty(idx), error('Quadrant errors of 100%, local bias is NaN'); end
            precPM=circstd(m(idx,8)-m(idx,6));
            varargout{1}=precPM;
            meta.ylabel='Local central precision error (deg)';
        case 'precPmedian'
            idx=find(abs(m(:,7))<=30);
            m=m(idx,:);
            precPM=circstd(m(:,8)-m(:,6));
            varargout{1}=precPM;
            meta.ylabel='Central precision error (deg)';
        case 'rmsPmedianlocal'
            idx=find(abs(m(:,7))<=30);
            m=m(idx,:);
            idx=find(abs(mynpi2pi(m(:,8)-m(:,6)))<90);
            rmsPlocal=sqrt(sum((mynpi2pi(m(idx,8)-m(idx,6))).^2)/size(idx,1));
            varargout{1}=rmsPlocal;
            meta.ylabel='Local central RMS error (deg)';
            % chance performance (FIXME: not working with exp_baumgartner2014!)
            %       tang = m(:,6);
            %       rang = min(tang):max(tang); % assume listener knows target range
            %       p = ones(length(rang),length(tang)); % equally distributed response probabilities
            %       chance = baumgartner2014_pmv2ppp(p,tang,rang,'PE');
            %       par.chance = chance;
            %       meta.chance = 'Chance performance';
        case 'rmsPmedian'
            idx=find(abs(m(:,7))<=30);
            m=m(idx,:);
            rmsP=sqrt(sum((mynpi2pi(m(:,8)-m(:,6))).^2)/size(m,1));
            varargout{1}=rmsP;
            meta.ylabel='Central RMS error (deg)';
        case 'SCC'
            trgt = m(:,[1 2]);
            resp = m(:,[3 4]);
            T = pol2dcos(trgt);
            R = pol2dcos(resp);
            sxy = T' * R;
            sxx = T' * T;
            syy = R' * R;
            scc = det(sxy) / sqrt(det(sxx) * det(syy));
            varargout{1}=scc;
            meta.ylabel='Spatial correlation coefficient';

        case 'gainLstats'

            [l.b,l.bint,l.r,l.rint,l.stats]=regress(m(:,7),[ones(length(m),1) m(:,5)]);
            varargout{1}=l;
            meta.ylabel='Lateral gain';

        case 'gainL'

            l = localizationerror(m,'gainLstats');
            varargout{1}=l.b(2);
            meta.ylabel='Lateral gain';

        case 'precLregress'

            l = localizationerror(m,'gainLstats');
            x = -90:90;
            y = l.b(2)*m(:,5) + l.b(1); % linear regression
            dev = m(:,7) - y; % deviation
            dev = dev(abs(dev)<=45); % include only quasi-veridical responses
            prec = rms(dev);
            varargout{1} = prec;
            meta.ylabel = 'Lateral scatter (deg)';

        case 'pVeridicalL'
            tol = 45; % tolerance of lateral angle error for quasi-veridical responses
            l = localizationerror(m,'gainLstats');
            x = -90:90;
            y = l.b(2)*m(:,5) + l.b(1); % linear regression
            dev = m(:,7) - y; % deviation wrapped to +-180 deg
            prop = sum(abs(dev)<=tol)/length(dev);
            varargout{1} = prop*100;
            meta.ylabel = '% lateral quasi-veridical';

        case {'sirpMacpherson2000', 'sirp'}
            delta = 40; % outlier tolerance in degrees
            Nmin = 5; % minimum number of responses used for regression
            maxiter = 100; % max number of iterations

            % Only responses for targets located within 30 degrees of the median plane are included
            m = m( abs(m(:,5))<=30,: );

            for frontback = 1:2 % Analyze front and back hemispheres separately

                if frontback == 1 % Front
                    mhemi = m( m(:,6)<=90 ,:); % frontal target locations
                    idinit = mhemi(:,8)<=90; % initialization with correct responses
                else % Back
                    mhemi = m( m(:,6)>=90 ,:); % rear targets
                    idinit = mhemi(:,8)>=90;  % initialization with correct responses
                end

                if length(idinit)<Nmin
                    r.b = nan(1,2);
                    r.stats = nan(1,4);
                else
                    iditer = false(length(idinit),maxiter);
                    y = mhemi(:,8); % response
                    x = mhemi(:,6); % target
                    for iter = 1:maxiter
                        B=regress(y(idinit),[ones(sum(idinit),1) mhemi(idinit,6)]);
                        yhat=B(2)*x+B(1);
                        dev = mod( (y-yhat) +180,360) -180; % deviation wrapped to +-180 deg
                        iditer(:,iter)= abs(dev) < delta;
                        if isequal(iditer(:,iter),idinit)
                            break % because it converged
                        else
                            idinit = iditer(:,iter);
                        end
                    end
                    [nselect,iopt] = max(sum(iditer));
                    idselect = iditer(:,iopt);

                    [r.b,r.bint,r.r,r.rint,r.stats]=regress(...
                        mhemi(idselect,8),[ones(nselect,1) mhemi(idselect,6)]);

                end
                varargout{frontback}=r;
            end

            par='SIRP results';

        case 'gainPfront'

            f = localizationerror(m,'sirpMacpherson2000');
            varargout{1}=f.b(2);
            meta.ylabel='Frontal polar gain';

        case 'gainPrear'

            [tmp,r] = localizationerror(m,'sirpMacpherson2000');
            varargout{1}=r.b(2);
            meta.ylabel='Rear polar gain';

        case 'gainP'

            [f,r] = localizationerror(m,'sirpMacpherson2000');
            varargout{1}=(f.b(2)+r.b(2))/2;
            meta.ylabel='Polar gain';

        case 'slopePfront'

            g = localizationerror(m,'gainPfront');
            varargout{1} = rad2deg(acos(1./sqrt(g.^2+1)));
            meta.ylabel='Frontal polar regression slope';

        case 'slopePrear'

            g = localizationerror(m,'gainPrear');
            varargout{1} = rad2deg(acos(1./sqrt(g.^2+1)));
            meta.ylabel='Rear polar regression slope';

        case 'slopeP'

            g = localizationerror(m,'gainP');
            varargout{1} = rad2deg(acos(1./sqrt(g.^2+1)));
            meta.ylabel='Polar regression slope';

        case 'pVeridicalPfront'

            tol = 45; % tolerance of polar angle error for quasi-veridical responses

            f = localizationerror(m,'sirpMacpherson2000');

            mf = m( abs(m(:,5))<=30 & m(:,6)< 90 ,:); % frontal central data only

            x = -90:270;
            yf=f.b(2)*mf(:,6)+f.b(1); % linear prediction for flat, frontal targets

            devf = mod( (mf(:,8)-yf) +180,360) -180; % deviation wrapped to +-180 deg

            pVer = sum(abs(devf) <= tol) / length(devf);

            varargout{1}=pVer*100;
            meta.ylabel = '% quasi-veridical (frontal)';

        case 'pVeridicalPrear'

            tol = 45; % tolerance of polar angle error for quasi-veridical responses

            [~,r] = localizationerror(m,'sirpMacpherson2000');

            mr = m( abs(m(:,5))<=30 & m(:,6)>=90 ,:); % rear central data only

            x = -90:270;
            yr = r.b(2)*mr(:,6) + r.b(1); % linear prediction for flat, rear targets

            devr = mod( (mr(:,8)-yr) +180,360) -180; % deviation wrapped to +-180 deg

            pVer = sum(abs(devr) <= tol) / length(devr);

            varargout{1}=pVer*100;
            meta.ylabel = '% quasi-veridical (rear)';

        case 'pVeridicalP'

            tol = 45; % tolerance of polar angle error for quasi-veridical responses
            [f,r] = localizationerror(m,'sirpMacpherson2000');

            mf = m( abs(m(:,5))<=30 & m(:,6)< 90 ,:); % frontal central targets
            mr = m( abs(m(:,5))<=30 & m(:,6)>=90 ,:); % rear central targets

            x = -90:270;
            yf = f.b(2)*mf(:,6) + f.b(1); % linear prediction for flat, frontal targets
            yr = r.b(2)*mr(:,6) + r.b(1); % linear prediction for flat, rear targets

            devf = mod( (mf(:,8)-yf) +180,360) -180; % deviation wrapped to +-180 deg
            devr = mod( (mr(:,8)-yr) +180,360) -180; % deviation wrapped to +-180 deg

            pVerf = sum(abs(devf) <= tol) / length(devf);
            pVerr = sum(abs(devr) <= tol) / length(devr);

            pVer = (pVerf + pVerr)/2;

            varargout{1}=pVer*100;
            meta.ylabel = '% quasi-veridical';

        case 'precPregressFront'

            tol = 45; % tolerance of polar angle error for quasi-veridical responses

            f = localizationerror(m,'sirpMacpherson2000');

            mf = m( abs(m(:,5))<=30 & m(:,6)< 90 ,:); % frontal central targets only

            x = -90:270;
            yf = f.b(2)*mf(:,6) + f.b(1); % linear prediction for flat, frontal targets

            devf = mod( (mf(:,8)-yf) +180,360) -180; % deviation wrapped to +-180 deg

            dev = devf(abs(devf)<=tol); % include only quasi-veridical responses
            prec = rms(dev);
            varargout{1} = prec;
            meta.ylabel = 'Frontal polar scatter (deg)';

        case 'precPregressRear'

            tol = 45; % tolerance of polar angle error for quasi-veridical responses

            [~,r] = localizationerror(m,'sirpMacpherson2000');

            mr = m( abs(m(:,5))<=30 & m(:,6)>=90 ,:); % rear central targets only

            x = -90:270;
            yr = r.b(2)*mr(:,6) + r.b(1); % linear prediction for flat, rear targets

            devr = mod( (mr(:,8)-yr) +180,360) -180; % deviation wrapped to +-180 deg

            dev = devr(abs(devr)<=tol); % include only quasi-veridical responses
            prec = rms(dev);
            varargout{1} = prec;
            meta.ylabel = 'Rear polar scatter (deg)';

        case 'precPregress'

            precF = localizationerror(m,'precPregressFront');
            precR = localizationerror(m,'precPregressRear');
            prec = (precF+precR)/2;
            varargout{1} = prec;
            meta.ylabel = 'Polar scatter (deg)';

        case 'perMacpherson2003'

            if isempty(kv.f) || isempty(kv.r)
                amt_disp('Regression coefficients missing! Input is interpreted as baseline data!','documentation');
                amt_disp('For this analysis the results from `sirpMacpherson2000()` for baseline data','documentation');
                amt_disp('are required and must be handled as `localizationerror(m,f,r,...)`.','documentation');
                [kv.f,kv.r] = localizationerror(m,'sirpMacpherson2000');
            end

            plotflag = false;
            tol = 45; % tolerance of polar angle to be not counted as polar error

            mf = m( abs(m(:,5))<=30 & m(:,6)< 90 ,:); % frontal central targets only
            mr = m( abs(m(:,5))<=30 & m(:,6)>=90 ,:); % rear central targets only

            x = -90:270;
            yf=kv.f.b(2)*mf(:,6)+kv.f.b(1); % linear prediction for flat, frontal targets
            yr=kv.r.b(2)*mr(:,6)+kv.r.b(1); % linear prediction for flat, rear targets

            devf = mod( (mf(:,8)-yf) +180,360) -180; % deviation wrapped to +-180 deg
            devr = mod( (mr(:,8)-yr) +180,360) -180;

            f.pe = sum(abs(devf) > tol);
            r.pe = sum(abs(devr) > tol);

            per = (f.pe + r.pe) / size(m,1) * 100; % polar error rate
            varargout{1} = per;
            meta.ylabel = 'Polar error rate (%)';

            if plotflag
                plot(m(:,6),m(:,8),'ko');
                axis equal
                axis tight
                hold on
                plot(mf(:,6),yf,'LineWidth',2);
                plot(mr(:,6),yr,'LineWidth',2);
                plot(mf(:,6),yf+45,'r','LineWidth',2);
                plot(mf(:,6),yf-45,'r','LineWidth',2);
                plot(mr(:,6),yr+45,'r','LineWidth',2);
                plot(mr(:,6),yr-45,'r','LineWidth',2);
                title(['PER: ' num2str(per,2) '%'])
            end

        otherwise
            error(['Unknown ERR: ' err]);
    end
    if length(varargout) < 2
        varargout{2}=meta;
    end
    if length(varargout) < 3
        varargout{3}=par;
    end
end

end

function out_deg=mynpi2pi(ang_deg)
ang=ang_deg/180*pi;
out_rad=sign(ang).*((abs(ang)/pi)-2*ceil(((abs(ang)/pi)-1)/2))*pi;
out_deg=out_rad*180/pi;
end

function [phi,mag]=circmean(azi,dim)
% [PHI,MAG]=circmean(AZI,DIM)
%
% Calculate the average angle of AZI according
% to the circular statistics (Batschelet, 1981).
%
% For vectors, CIRCMEAN(AZI) is the mean value of the elements in AZI. For
% matrices, CIRCMEAN(AZI) is a row vector containing the mean value of
% each column.  For N-D arrays, CIRCMEAN(azi) is the mean value of the
% elements along the first non-singleton dimension of AZI.
%
% circmean(AZI,DIM) takes the mean along the dimension DIM of AZI.
%
% PHI is the average angle in deg, 0..360�
% MAG is the magnitude showing the significance of PHI (kind of dispersion)
%   MAG of zero shows that there is no average of AZI.
%   MAG of one shows a highest possible significance
%   MAG can be represented in terms of standard deviation in deg by using
%       rad2deg(acos(MAG)*2)/2
%
% [PHI,MAG]=circmean(AZI,W) calculates a weighted average where each sample
% in AZI is weighted by the corresponding entry in W.
%
% Piotr Majdak, 6.3.2007

if isempty(azi) % no values
    phi=NaN;
    mag=NaN;
    return;
end

if prod(size(azi))==1 % one value only
    phi=azi;
    mag=1;
    return;
end

if ~exist('dim')  % find the dimension
    dim = min(find(size(azi)~=1));
    w=ones(size(azi));
elseif prod(size(dim))~=1
    % dim vector represents weights
    w=dim;
    dim = min(find(size(azi)~=1));
end


x=deg2rad(azi);

s=mean(sin(x).*w,dim)/mean(w);
c=mean(cos(x).*w,dim)/mean(w);

y=zeros(size(s));
for ii=1:length(s)
    if s(ii)>0 & c(ii)>0
        y(ii)=atan(s(ii)/c(ii));
    elseif c(ii)<0
        y(ii)=atan(s(ii)/c(ii))+pi;
    elseif s(ii)<0 & c(ii)>0
        y(ii)=atan(s(ii)/c(ii))+2*pi;
    elseif s(ii)==0
        y(ii)=0;
    else
        y(ii)=NaN;
    end
end
mag=sqrt(s.*s+c.*c);

phi = mod(rad2deg(y), 360);
phi((phi==0) & (rad2deg(y)>0)) = 360;

end

function mag=circstd(azi,dim)

% MAG=circstd(AZI,DIM)
%
% Calculate the standard deviation angle (in deg) of AZI (in deg) according
% to the circular statistics (Batschelet, 1981).
%
% For vectors, CIRCSTD(AZI) is the std of the elements in AZI. For
% matrices, CIRCSTD(AZI) is a row vector containing the std of
% each column.  For N-D arrays, CIRCSTD(azi) is the std value of the
% elements along the first non-singleton dimension of AZI.
%
% CIRCSTD(AZI,DIM) takes the std along the dimension DIM of AZI.
%
%
% MAG=CIRCSTD(AZI,W) calculates a weighted std where each vector strength
% of AZI is weighted by the corresponding entry in W.
%
% Piotr Majdak, 21.08.2008

if isempty(azi) % no values
    mag=NaN;
    return;
end

if prod(size(azi))==1 % one value only
    mag=0;
    return;
end

if ~exist('dim')  % find the dimension
    [phi,mag]=circmean(azi);
elseif prod(size(dim))~=1
    [phi,mag]=circmean(azi,dim);
end

mag=rad2deg(sqrt(2*(1-mag)));

end